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Forward 

Driving simulation requires the construction of a total synthetic real¬ 
time environment. An effective simulator must provide proper mo¬ 
tion, visual, sound, and tactile cueing for the driver to respond in 
a natural manner. A simulated driving environment is an excellent 
example of a complex “real-time microworld”, where, a person must 
interact completely and naturally with an artificial, but compelling, 
computer controlled world. 

Driving simulation is an interesting case study because it entails 
synthesizing a familiar, yet highly non-trivial, environment. It must 
provide accurate vehicle dynamics predicting real vehicle behavior; 
graphics which provide sufficient optical flow density for proper con¬ 
trol and speed-distance judgments by the subject; optimal accelera¬ 
tion cueing; sound and control force feedback cueing mechanisms to 
enhance the sense of realism for the subject. The goal of human inter¬ 
action is to provide a designer with information about performance, 
stability, ride quality, and safety margins to suit human tolerances in 
control, task loading, and comfort. 

Our purpose in writing these notes is to highlight the essential sys¬ 
tem issues for the design and implementation of real-time, animated 
environments for vehicle (and, in particular, driving) simulation. 

To meet these goals, these notes introduce ideas about 

• Total simulator design. 

• Real-time vehicle dynamics. 

• Host computer and parallel processing issues. 

• Physiological models of perception for cueing optimization. 

• Washout algorithm design for matching motion system capabil¬ 
ities. 

• Control force loading. 

To know who to blame, the first author (R.D.) was primarily re¬ 
sponsible for the first three chapters on dynamics and the real-time 



6 


vehicle algorithm, while the second author (D.I.) was primarily re¬ 
sponsible for chapter four on physiological considerations. 

We want to thank the members of the Driving Simulation Project 
at Evans & Sutherland for their efforts towards real-time vehicle sim¬ 
ulation. Pete Doenges and John Briggs have provided the vision, 
conceptual foundations for advanced driving simulation, and part of 
the text from their Siggraph paper ([22]). We also want to thank E. 
Haug, D-S. Bae, R-S. Hwang, D. Andersen, and W-S. Lee. Any errors 
in these notes are, of course, our own contribution. 



Chapter 1 

The Driving Simulator 
System 


1.1 Introduction 

In the last fifteen or so years, man-in-the-loop simulation using computer- 
generated, real-time graphics has assumed a major role not only in 
pilot training, but also in engineering applications for land, sea, air, 
and space vehicles. Man-in-the-loop simulation uses a person in the 
control loop of a simulation to provide feedback to the system opera¬ 
tions. The operator experiences a synthetic environment, manipulates 
his controls in response to the situation, and experiences the results 
of his actions. In a real-time simulator, the dynamics model is con¬ 
tinually evaluated to give feedback to the operator. On the basis of 
this data, the simulator provides realistic environmental cues to allow 
him to control the simulated vehicle in the same way that he would 
in the real world. 

The goal of human interaction with a design is to tune raw per¬ 
formance, stability, ride quality, and safety margins, to suit human 
tolerances in control, task loading, and comfort. Through the use of 
simulation, training or design testing can often be achieved at lower 
cost, experiments can be conducted under controlled, repeatable cir¬ 
cumstances, and — important for engineering — vehicle designs can 
be evaluated while minimizing the number of prototypes. 

An important engineering use of simulation with a person in the 
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loop is to determine vehicle response and to study the man-machine 
interface before the start of a vehicle test program. For example, a 
new aircraft can be modeled and test flown even before it leaves the 
drawing board. The pilots can give valuable feedback to the designers 
about the plane’s characteristics at an early enough stage to permit 
changes in the airframe. Integration of the avionics systems can be 
achieved at an earlier time. Workload measurements can be made, 
and different alternatives can be tested. Computer graphics represents 
a large portion of the cost of man-in-the-loop simulation. Hence, as 
the capabilities of computer image generation have improved and the 
relative cost of simulation has dropped, there are better economic 
justifications for its use. In some applications, such as the space shut¬ 
tle simulator, there has never been a good alternative to simulator 
training. In other areas, however, such as vehicle engineering, cost 
of man-in-the-loop simulation testing has always been an important 
consideration, and, as a result, has only recently been low enough to 
justify simulation on a large scale. 

Computer graphics is the critical element in creating a realistic 
operator environment in a simulator, but the graphics must work in 
concert with the rest of the system to be believable. In particular, the 
vehicle motion and the scenes generated must seem realistic in order 
not to seem like an arcade game. To show how real-time computer 
graphics and dynamics must work together in vehicle simulation, we 
have selected a specific application to discuss — driving simulation. 
We will give a brief history of the use of synthetic imagery in simu¬ 
lation, identify the components of a modern simulator with emphasis 
on the role of the visual system, present a new approach to computing 
vehicle dynamics that promises significantly improved vehicle models 
for realistic motion, and conclude with a view of the future. 


1.2 Historical Perspective 

Department of Defense initiatives through the military and its con¬ 
tractors have done much to advance the state-of-the-art in simula¬ 
tion and computer graphics. For example, seminal work in computer 
graphics conducted at the University of Utah under the leadership 




of Dr. David Evans two decades ago was largely sponsored by the 
Department of Defense’s Advanced Research Projects Agency. 

The U. S. Air Force supported in part the first use of CRTs in 
driving simulation at MIT in 1962 [63], This work by T.B. Sheridan, 
H.M. Paynter, and S.A. Coons was the earliest use of an artificial 
display of a vehicle path in dynamic perspective. A refined imple¬ 
mentation was developed by W.W. Wierwille, G. A. Gagne, and J. 
Knight at Cornell Aeronautical Laboratory in 1966 (69]. 

Use of these same techniques occurred at General Motors in the 
late 1960’s. In the early 1970’s, both Volkswagen [59] and the Swedish 
Road and Traffic Institute (VTI) [52] built driving simulators with vi¬ 
sual and motion subsystems. A number of simulators were built over 
the next decade. For example, simulators at the Virginia Polytechnic 
Institute and State University [17] and Deere and Company [26] were 
developed primarily to support human factors research. Another sys¬ 
tem, complete with motion platform, was constructed at Universitaet 
der Bundeswehr Muenchen for research into guiding land vehicles with 
computer vision — a sort of robot-in-the-loop simulation [23]. The 
graphics associated with these systems were limited. Line segments 
representing the road were displayed on a single CRT or projected on 
a flat screen in front of the driver. 

In 1984, over twenty years after the MIT project, Daimler-Benz 
built, in Berlin, the first driving simulator expressly as an aid to 
designing vehicles [27] [33]. This system marked the first use of a 
180-degree-wide image on a dome in a driving simulator. For the first 
time a driver was able to use his peripheral vision in a natural way. 
The Daimler-Benz simulator was built to allow test drivers to make 
subjective judgments about the handling qualities of new automobile 
designs. 


1.3 Today’s Advanced Simulator 

An advanced man-in-the-loop simulator vehicle (Figure 1.1) is made 
up of a few, expensive subsystems that combine to create a complete 
synthetic driving environment for the operator: 



10 


CHAPTER 1. THE DRIVING SIMULATOR SYSTEM 



Figure 1.1: Modern advanced vehicle simulator system 
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• The host computer (Figure lg) controls all of the real-time op¬ 
erations of the system, including computation of the vehicle and 
motion platform dynamics. All calculations for the entire sys¬ 
tem are performed within a 10 ms time frame. An integrated 
data base contains visual models, geometric data on the simu¬ 
lated environment, data relating to the sounds encountered, and 
mechanical characteristics of the vehicle. 

• The visual system calculates the correct image (Figure lb) based 
on position information from the host computer. It removes 
hidden areas, computes the sun angle, smooth shades, textures, 
and anti-aliases the result. There is a high degree of parallelism 
in the foregoing operations, allowing about ten GIPS (Giga In¬ 
structions Per Second) operation in a large system. 

• Multiple CRT projectors (Figure la) display an image on the 
dome. Special screen material is used for maximum brightness 
and contrast. 

• The motion platform (Figure 1-2,3,d) is a configuration of hy¬ 
draulic actuators that can move the operator’s cab in any direc¬ 
tion or angle by combining motions about the x, y, z, yaw, pitch, 
and roll axes. The control signals to the motion platform are fil¬ 
tered (“washed out”) to transform the virtual real-world motion 
computed by the dynamics software into movements achievable 
by the system. 

• The sound subsystem (Figure le) uses multiple speakers to gen¬ 
erate appropriate sounds to simulate the engine, gear whine, 
wind noise, passing traffic, etc. 

• Instrumentation (Figure lh) provides data to the host computer 
to monitor the the entire system and ensure the driver’s safety. 

• Control force loaders (Figure lc) are used to provide appropriate 
forces back to the operator through the steering wheel, brake 
pedal, and transmission lever to replicate the feedback forces 
found in a real car. 
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• The operations console (Figure If) allows a single person to 
oversee all of the simulation activities, allowing flexible control 
over vehicle tests. 


1.4 Relating Graphics and Dynamics Un¬ 
der Real-Time Constraints 

Driving simulation for automotive engineering must provide enough 
accuracy in simulating vehicle motion and in subsystems cueing (vi¬ 
sual, sound, motion, control force feedback, and instruments) to stim¬ 
ulate realistic driver responses that accurately portray flaws or inad¬ 
equacies in the “paper” design of the target vehicle. The fidelities of 
graphics and dynamics crucial to driver cueing are measured by dif¬ 
ferent yardsticks than those employed in non-real-time modeling. For 
example, the graphics realism of ray-tracing or total correlation be¬ 
tween dynamics response from instrumented cars and vehicle models 
needing overnight computer runs are not necessarily pertinent. 

Fully interactive real-time simulation demands care in the relative 
realism of the environmental scene (scene specificity and density) and 
dynamics accuracy significant to the driver response (chassis acceler¬ 
ations and control forces). Also important is the rate at which new 
cues are delivered to the driver and the time lag between his control 
commands and the corresponding cue changes. What is the proper 
balance between real-time constraints and the cueing rates needed to 
support a driver’s decision process? What is the allowable system 
time lag before false driver responses are induced? 

The subtlety of human perception and response while interacting 
with a vehicle forces difficult compromises. The dynamics consid¬ 
erations bearing on driver behavior demand an almost open-ended 
fidelity escalation. A driver’s sensitivities to visual and tactile cues 
tend to peak in the region of a few hertz while steering in response to 
sudden external stimuli and at tens of hertz for periodic and impulsive 
road features. Rates of cueing tend to cap the response bandwidth of 
simulator vehicle dynamics. Much higher frequencies can arise from 
mechanical resonances of suspension parts and body torsion/flex, or 
the dynamic response of an on-board digital control system. How will 
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the dynamics simulation deal with cueing-significant modal energy 
which is coupled by system non-linearities from these higher frequency 
ranges down into the ranges of driver sensitivities? 

Simulation-aided engineering requires vehicle designers to trust 
simulation realism when no equivalent physical vehicle system has 
been prototyped or tested. Computer-aided design and analysis tools 
can often supply reliable subsystem data. Since the total system can 
behave differently from the sum of its parts, how can the designer 
trust a simulator whose response might compromise the otherwise 
adequate subsystem data? 

The answers to all these tough questions can be summarized with 
a few rules of thumb: 

• Provide flexibility in trading system resources to enhance tem¬ 
porarily the specific areas of experimental focus at the expense 
of others. 

• Ensure that the minimum requirements of human perceptual 
cueing thresholds and of vehicle operations are met by simulator 
performance and functionality. 

• Attempt to drive overall simulator response as close to the step¬ 
steering and steady state (in both phase and frequency) re¬ 
sponses of the vehicle as practical. 

• Provide constant interfaces between different forms of the same 
data distributed within the simulator so that changes afforded 
by technology advances can be easily incorporated into simula¬ 
tor improvements. 

• Develop simulator software architectures that port easily into 
future computer systems. 

• Provide as much reserve capability as economically possible. 
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1.5 CIG Demands of High-Performance 
Driving Simulation 

The use of real-time, man-in-the-loop simulators for the evaluation 
of designs of high-performance land vehicles has greatly lagged flight 
simulation. Major advances in dynamics, washout, motion, and graph¬ 
ics have recently set the stage for another leap. 

The current state of CIG (Computer Image Generation) evolution 
is closely approaching photographic realism with update rates, lag, 
and image quality that suffice for a broad range of operational and 
perceptual objectives. Recent advances in computing and mechanical 
dynamics techniques have set the stage to achieve vehicle behavior in 
a simulator rivaling the quality of CIG real-time graphics, and to af¬ 
ford dynamics quality capable of unmasking subtle CIG shortcomings 
rarely encountered before now. 

CIG improvement is driven by operational maneuvering and per¬ 
ceptual requirements (24). Driving involves very low eye heights where 
optical flow density must change very rapidly over the field of vision 
available in a car. The driver must be able to judge speed and prox¬ 
imity to obstacles very quickly by visualizing textural cues in and 
around the road as well as passing 3D features. The visual depiction 
of the road surface must be good enough for the driver to anticipate 
road hazards like bumps and ice. When appropriate, urban, coun¬ 
try, and highway scenes must be possible with contending traffic and 
pedestrians to saturate the driver’s task loading and reaction times. 

The basic psychology of vision and kinematics of driver responses 
at the controls of the car also demand CIG fidelity. The effective 
resolution of the CIG image determines what the driver can detect 
and recognize. A driving simulator should ideally match the lag, from 
control input to visual result, at a few tens of milliseconds typical of 
responsive cars. The rate of image refresh and update must satisfy 
eye demands for visual tracking, smooth motion of scenery, flicker-free 
viewing and discrimination of temporal events that may alter driver 
response. 

The current generation of top-end real-time CIG is represented by 
such products as the Evans & Sutherland ESIG-1000. This system 
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Figure 1.2: Advanced real-time display. Example of high level of 
detail in modern real-time image generation using textures. 

olTers critical improvements in scene performance needed to discrimi¬ 
nate the nuances of vehicle behavior possible with the dynamics tech¬ 
niques covered in the following section. Synthetic and photographic 
texture along with high 3D feature densities now allow sound judg¬ 
ments to be made about speed and distance. Subpixel anti-aliasing 
allows road signs and road surfaces to be visualized reliably at varied 
vehicle speeds. Generalized occultation allows complex traffic move¬ 
ment over rolling terrain, around urban clutter, and among other 
traffic. Multiple eyepoints can support the virtual eyepoints of mir¬ 
rors and the simultaneous observation of own-vehicle chassis and sus¬ 
pension behavior from external vantage points while driving the car 
through difficult maneuvers. 

Figure 1.2 shows a typical driving landscape and situation achiev¬ 
able with ESIG-1000 for driving simulation. Figure 1.3 represents 
the kind of car-external viewing of critical moments in suspension de¬ 
formation and chassis excursions. This real-time visual display can 
complement a designer’s study of subtle car subsystem motion over 
specific road geometries. 
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Figure 1.3: Examples of real-time generation of an external view of 
the vehicle motion showing suspension deformation. Frames are 80 
ms apart. 




Chapter 2 
Dynamics 


2.1 Dynamics for Computer Simulation 

This chapter of the the notes shows how to make complex multi-body 
vehicles in the microworld move correctly in real-time. It reviews the 
basic concepts of rigid body kinematics and dynamics, including the 
virtual work equations and Lagrange multipliers. The next chapter 
will apply these ideas in the form of Bae’s recursive algorithm for 
linked systems of rigid bodies. This algorithm has been shown to 
be the fastest available for either serial or shared-memory parallel 
processors, and its description is the goal of the dynamics part of 
these notes. 

The following topics are covered in the next two chapters dealing 
with real-time vehicle dynamics and its implementation for a driving 
simulator: 

• Quick review of basic particle and rigid body dynamics. 

— Newton’s Laws. Linear and angular momentum. 

— Virtual work principle. Constraints. Generalized coordi¬ 
nates. 

- Rigid body kinematics. Transformation matrices. Euler 
angles and parameters. 

— Newton-Euler equations for a rigid body. Example of sym¬ 
metric top. Constrained virtual work equations for systems 


17 



18 


CHAPTER 2. DYNAMICS 


of rigid bodies. 

• Real-time dynamics algorithms. 

— Use of system graph to decompose linked rigid body vehi¬ 
cles. Example of car. 

— Description of Bae’s recursive algorithm for real-time dy¬ 
namics of linked, multi-body systems. Example of car. 

— Parallel processing issues. Alliant FX/80 example. 

• Practicalities. It will discuss the following topics: 

— Real-time integration methods. Using an eigenvalue anal¬ 
ysis. 

- Force models. 

Newton started by considering the motion of planets and projec¬ 
tiles. Later investigators considered rigid bodies. Modern dynamics 
is mainly concerned with complex systems of linked, constrained rigid 
and flexible bodies, moving under the influence of external forces and 
moments. Constrained multi-body dynamics is an active area of cur¬ 
rent research. Its roots go well back to to the late eighteenth and early 
nineteenth centuries. The theoretical foundations are relatively well 
understood, with many possible approaches. For complex systems, 
few analytic models could (or can) be solved. With the advent of 
computers, involved models can now be done easily in non-real-time. 
The current challenge is to solve these same models in real-time to 
allow man-in-the-loop interactions. 


2.2 Review of Newtonian Dynamics 

In this section we review the concepts from classical mechanics needed 
to understand the vehicle dynamics algorithms. 

2.2.1 Newton’s Laws of Motion 

Newton’s laws originated in considering planetary motion. An indi¬ 
vidual planet is a body moving in space subject to various known 
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gravitational forces from the sun, other planets, and any other bodies 
having mass. Newton’s original concern with multi-body dynamics 
was to derive Kepler’s laws of planetary motion from first principles, 
i.e. Newton’s law of universal gravitation. One of his first theorems 
showed that a planet was, for purposes at hand, equivalent to a par¬ 
ticle with the total planet mass. 

Our job is more complex, since we must deal with bodies which 
do not act like point particles; our bodies will have rotational, as well 
as translational, inertia about their center-of-mass. In addition, these 
bodies are linked in complex kinematic configurations, where relative 
motion is constrained by the presence of joints. We’ll defer until a 
later section a full treatment of the rotational inertia of a rigid body. 
For now, a rigid body can be treated dynamically as a point particle, 
with all its mass at its center-of-mass and no rotational degrees-of- 
freedom (cf. 2.3). 

We are assuming that the reader is generally familiar with New¬ 
ton’s three laws and their consequences. 

2.2.1.1 Coordinate Systems 

It is a general principle that the results of physical experiments should 
be independent of the coordinates with which they are expressed. In 
particular, we should be able to describe the motion of multiple bodies 
within a vehicle by any convenient choice of coordinates; not just by 
the global cartesian coordinates of each body. The first part of these 
notes makes use of three-dimensional cartesian vector coordinates r = 
(x, y, z) with respect to some fixed global inertial coordinate frame 
(cf. Section 2.2.1.3). This is shown in Figure 2.1. In later sections 
we will show how to expand our horizons, making use of the simplest 
“generalized” coordinates to express our equations of motion. 

2.2.1.2 Newton’s Second Law 

If a “particle” has mass m and position vector r in a suitable coor¬ 
dinate system, then the first derivative with respect to time of its 
position vector is its velocity v ( so that v — r) and its second deriva¬ 
tive is its acceleration a (thus a = F). The linear momentum of the 
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Figure 2.1: Cartesian coordinate system for a particle, 
particle is the vector p = mv. 

We now can present the only dynamics equation you will ever need: 
Newton’s Second Law states that 

’p = F (2.1) 

where F is the vector sum of all forces acting on the body. If the 
mass is constant with time (we assume this is always the case in these 
notes), this is usually stated as ma = F. 

Everything we do is now downhill from this lofty summit. In 
essence, all other dynamics formulation are mere restatements of New¬ 
ton’s second equation; where the attempt is made to make it more 
convenient. For example, as we will see later, often it is very difficult 
to know the reaction forces which are generated by constrained mo¬ 
tion. Newton’s equation in the form of Equation 2.1 is then replaced 
by the more convenient virtual work equations. 


Example of Newton’s Second Law: Vertical motion of a mass m under 
gravitational acceleration g. The second law gives m'z — —mg (note the use 
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in two different ways of Newton’s law - if we know the constant gravitational 
acceleration g, the gravitational force must be mg). Integrating this equation 
gives z(t) = —gt 2 /2 + i 0 t 4- z$ where the initial conditions for velocity and 
position provide zq and zq. If the mass has initial height z — 0 and velocity 
io = vo then z(t) = —gt 2 /2 + vot. Because Newton’s equation is second- 
order, we must always prescribe both an initial position and velocity in order 
to solve the equation. 


2.2.1.3 Inertial Reference Frames 

We must remember that Newton’s Second Law, Equation 2.1, is valid 
only in inertial reference frames. These are the special frames which 
are in uniform motion with respect to one-another. Only in these 
frames do the equations of dynamics become so simple. For our ve¬ 
hicle models we use an “inertial” Cartesian reference frame fixed to 
the earth 1 . The final equations of vehicle motion are most conve¬ 
niently written in terms of the earth’s “inertial” frame, and not with 
respect to, say, a frame fixed on the vehicle chassis. Such a frame 
would not remain, in general, inertial, and the equations of motion 
would have additional funny force terms proportional to mass due to 
accelerations. 

2.2.1.4 Newton’s First Law 

Notice that Newton’s equation 2.1 tells us that if F = 0, then we 
have a = 0. Since a = r , we see that r = constant, and the body’s 
momentum and velocity are constant. Thus we have Newton’s First 
Law to the effect that when forces are not acting, a body remains tn 
uniform motion moving in a straight line. The law of conservation of 
linear momentum states that in the absence of external forces, linear 
momentum is conserved (cf. the next paragraph). Linear momentum 
is then the tendency to continue moving in the same direction. 

1 An alert reader will see that such a frame is not inertial due to the earths 
rotation, but the effects are so small that to an engineering approximation we can 
treat it as inertial for a vehicle on a road. 
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2.2.1.5 Newton’s Third Law 

Newton gave yet another law to complete his dynamics. It is a 
statement about linear momentum conservation when internal forces 
acting along straight lines between bodies are present. Newton’s 
Third Law state that to every action, there is an equal and opposite 
reaction 2 . 


Example of Newton’s Third Law: A rocket moves forward because it 
ejects, due to chemical reactions generating heat, high-speed gases towards 
the rear. The rocket motion is the “reaction” to the “action” of the gas 
motion. 


An important application of this law is to systems of particles. A 
system of particles without external forces acting cannot change its 
total linear momentum vector P — YliPx- This allows us to eliminate 
the (possible) internal forces between each pair of particles, and show 
that the motion of the total system acts as if all the system’s mass 
was concentrated at a single point (the center of mass), which, when 
no external forces act, itself has a uniform motion. 


Example of conservation of linear momentum: A bomb with mass m 
moving forward at velocity v explodes, the total system momentum must 
remain the same, so P = p, remains constant even though individual 
bomb fragments have many different momenta. 


2.2.1.6 Angular Momentum 

—t 

We remember that the angular momentum L of a particle is defined 
as the “moment of liner momentum”: 

L = r x p (2.2) 

2 This fails to be true in for some forces in electrodynamics since the magnetic 
Lorentz force does not act in this way. 



Figure 2.2: Angular momentum of particle moving in a circle. 


where p is the particles linear momentum, and r is the position or 
radius vector from the fixed origin of coordinates to the particle. An¬ 
gular momentum is intuitively considered as the tendency to continue 
rotating about a given point. In general, the angular momentum will 
clearly depend upon the choice of the center of coordinates. The 
torque N about the center of coordinates is the “moment of force F” 
is defined to be 


N = r x F 


(2.3) 


Example of angular momentum: A particle with mass m moving in a 
circle of radius r with uniform velocity v has an angular momentum of L = 
rp = mrv about the center in the polar direction normal to the plane of 
motion. If w is its angular velocity, then v = ru> and L = mr 2 w. Figure 2.2 
shown this example. 
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2.2.1.7 Rate of Change of Angular Momentum 

If we now use Newton’s second law, equation 2.1, and take the x- 
product of both sides with the radius vector r, we get 

N = r x j) (2.4) 

Using the definition of the angular momentum vector L , this can be 
rewritten as 

N = L (2.5) 

We see that the rate of change of a particle’s angular momentum 
is equal to the applied torque. If no torques are present, then the 
particle’s angular momentum is constant. This is the law of conser¬ 
vation of angular momentum. In analogy with Newton’s Third 
Law, for a system of particles, the total angular momentum L = L,- 

is conserved in the absence of external torques. 


Example of conservation of angular momentum: For a central force, 
such a gravity, which acts on a line between the sun and a planet, the torque 
N = f x r = 0, and so the angular momentum of the planet , and not just 
the total sun-planet system, is conserved. 


2.2.2 Multiple Particles and Constraints 

As an introduction to multi-body dynamics, we first look at multi- 
particle systems. Since a set of free or interacting particles is not 
interesting for our ultimate application, we consider the existence of 
constraints restricting the free motion of a particle set. 

2.2.2.1 Holonomic Constraints 

A constraint is some extra condition between the particles’ coordi¬ 
nates that limits the motion of the particles in some way. 


Examples of constraints: 



Figure 2.3: Constrained system - a Macpherson front suspension. 


1. A bead on a wire is constrained to one dimensional motion by the 
wire. 

2. A tire is constrained to revolve around the axis of its revolute joint. 

3. A Macpherson front suspension is constrained to have only a vertical 
“jounce” and horizontal “steer” motion. The steer is itself constrained 
by the tie-rod to be a function of steering wheel position. This is shown 
in Figure 2.3. 


A holonomic constraint is a explicit (possibly time dependent) 
relation among the particle coordinates r, that can be expressed in 
the form 

$(ri,r 2 ,...,t) = 0 (2.6) 

Other types of non-holonomic constraints exist (see [30]), but the 
holonomic ones are the most important, as well as the most pleasant 
to deal with in practice, since they reduce the dimensionality of the 
problem. 
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Examples of holonomic constraints: Various mechanical joint such as 
translational, revolute, and spherical joints all give holonomic constraints 
since they can all be expressed in terms of coordinate conditions, without 
using differentials. In addition, distance constraints are also holonomic. 

1. Once a revolute joint axis n is specified, any motion of a vector $ 
revolving about the joint must preserve the angle 9 = n • £/||n||||*||. 
The holonomic constraint is then $ = n • */||«||||*|| — cos(0) = 0. 

2. The existence of a tie-rod between steering rack and Macpherson sus¬ 
pension spindle enforces a holonomic distance constraint between the 
tie-rod attachment vectors. 


Example of a non-holonomic constraint: In one dimension, if a tire rolls 
without slipping we have * = r0, so = dx — rd9 is an integrable relation. 
See Figure 2.4. Thus we have a simple holonomic constraint, and we need 
only use either x or 9 as the coordinate which describes the system. In two 
dimensions, the relation between position and angle for a tire rolling with¬ 
out slipping on a plane is not integrable. This constraint involves relations 
between the coordinate differentials. Let x,y be the coordinates of the cen¬ 
ter of the wheel, let 9 be the angle the wheel has rotated, and let <j> be the 
angle between the disk and the x axis. Then since the wheel rolls without 
slipping, the center velocity must be related to rotation rate by v = t9 and 
the components of center motion are 

x = t>sin(<£) 

y = -vcos(<£) 

These give the differential relations between the variables 

dx - rd<f>sin(<j>) = 0 
dy + rd«f> cos(<f>) = 0 

It is not possible to integrate this differential relation to get a holonomic 
constraint between tire hub coordinates and tire rotation angle. The non- 
holonomic constraint gives conditions between the differentials which must 
be satisfied, without reducing the number of coordinates required to describe 
the motion. 
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Figure 2.4: Constrained system - wheel rolling without slipping. 


2.2.2.2 Multiple Constraints 

Let us assume we have M particles each with coordinates rj-. Of¬ 
ten such particles are not free and subject only to known forces, but 
also move under the influence of reaction forces provided by the con¬ 
straints. The existence of reaction forces are what force the system to 
obey the k constraint equations given by some large set of holonomic 
constraints such as 


$*(?!, r 2 ,...,r M ,0 = 


0 

0 (2-7) 

0 


For notational purposes, it is easier to write the Cartesian coordinates 
of all the bodies as the components of one large vector x of length N, 
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I 

where 

Xl 

= f l.* 


Xi 

= 


x 3 

= r M 


X 4 

= r 2(X 


X N 

II - 

N 


( 2 . 8 ) 


Using this notation we can consider the above set as k components of 
a vector constraint equation given by 


$(x) = 0 


(2.9) 


2.2.2.3 Constraint Jacobian 

The set of holonomic constraint conditions determines some maximal 
amount of dependence among the coordinates r,-. The amount of 
dependence is measured by the row-rank of the constraint Jacobian 
matrix whose ij-th component is 


< 6 - = 



( 2 . 10 ) 


where 1 < i < k and 1 < j < N. The row-rank is the number of 
linearly independent rows of the matrix, and is thus the dimension of 
the image space for the associated linear transformation. Something 
special happens when the row rank is as large as possible; that is, 
when it is equal to the number k of constraint equations. This occurs 
when all the constraints are independent and thus by the implicit 
function theorem 3 , we can split our original coordinates x into a of 
N-k independent coordinates and a set of k dependent coordinates. 
Thus, after renumbering the coordinates, we have 

x l = fl ( x flf—ki •f'N — fc+1 > • ■ •) 

x 2 = /2(ZN-*,Z/V-H-1>*--) , , 

. I 2 - 11 ; 

Xk = fk{XN-ki x N-k+U • • 0 


3 See any good advanced calculus text or look in [37]. 
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This is how independent constraints can be used to eliminate coor¬ 
dinate variables. We will see this again when we discuss generalized 
coordinates. 

An effectively computable method of finding the independent coor¬ 
dinates is to take the LU-decomposition or, better yet, the singular 
value decomposition 4 of the Jacobian. 


Example of a constraint Jacobian of two linked particles: Consider two 
particles moving on the plane z = 0, which are also linked rigidly together 
at unit distance. We have Fj, F 2 as position vectors. Thus we have six 
coordinates Since the particles are the plane x$ = z\ = 0 and *6 = *2 = 0. 
Since distance between them is constant, we have (Fj - F 2) 2 = 1. Then the 
constraint Jacobian is given by 

( 0 0 1 0 0 0 \ 

0 0 0 0 0 1 (2.12) 

2xi 2x2 0 -2x4 -2x5 0 J 

which clearly has the maximal row rank of three since both particles cannot 
simultaneously be at the origin. 


In order to see why the Jacobian enters, we take equation 2.9 and 
form its differential for each component i: 

£ = 0 ( 2 - 13 ) 

j VX) 

This is just a set of linear equations in the dxf s with coefficients 
given by the partials derivatives of in other words, by the Jaco¬ 
bian matrix. The maximal set of independent differentials, and so of 
coordinates, is then given by the row rank of the Jacobian matrix. 

2.2.3 Principle of Virtual Work 

In this section, the subscript t with denote an individual particle with 
position vector r,. When the position vectors of a system of particles is 


4 See the reference (29|. 
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considered, we write r = (rj, r 2 ,...). This convention will be followed 
for the remainder of the section. 

3.2.3.1 Statement of Virtual Work Principle 

The Principle of Virtual Work states that for a dynamical system 
to be in equilibrium, the sum of the internal and external virtual work 
must vanish. Clearly, this deserves some elaboration. What is being 
said is that a mechanical system is in its equilibrium state if when, 
holding time fixed, we vary any of the coordinates infinitesimally in 
a way compatible with all constraints, no work is performed. For 
static equilibrium, only the external forces have to be included in 
the virtual work since the reaction forces from the constraints are in 
general workless. If dynamic equilibrium is wanted, we must include 
the work due to inertial forces. 

2.2.3.2 Relation with other Variational Methods 

We remark that the virtual work formulation of dynamics is the most 
general available. When only conservative forces are present 5 , then for 
static equilibrium, it is equivalent to the minimization of the system’s 
potential energy E = <j> 6 But no such potentials exist in general for 
dissipative forces such as viscous friction. An even more fundamental 
problem exists if the forces of constraint can themselves do work - as 
in the case of rolling friction. The total system energy is now total 
energy = internal energy - external work. If this is minimized, we 
recover our original virtual work principle. 

In general, for complex mechanical systems with dissipative el¬ 
ements, the virtual work principle seems in practice to be superior 

5 If a force can be written as the gradient of a potential fnnction f = V<f>, work 
done by the force does not depend on the path taken in configuration space, only 
the starting and ending positions matter. In this case the force is called called 
conservative. The work done by a non-conservative force, such as friction, clearly 
does depend upon the path. If we move a particle against friction and then return 
to the starting point, work has been done. 

°The kinetic energy of a particle is T = p 2 /2m, and it measures the energy of 
motion. Potential energy is a function only of position, e.g. height in a gravitational 
field. Only the sum of these two energies is conserved when no work is performed 
on or by the system. 
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to other less general variational techniques. Powerful methods due 
to Lagrange and Hamilton unfortunately are currently unsuitable for 
vehicle simulation. 


2.2.3.3 Derivation of Virtual Work Principle 

A virtual (infinitesimal) displacement of a system refers to a change 
in its configuration as the result of any arbitrary infinitesimal change 
of coordinates £r t , consistent with the forces and constraints imposed 
on the system at time t. The term virtual is used to distinguish these 
displacements form real kinematic displacements during the time dt 
when the forces and constraints may be changing. 

Suppose that the system is in equilibrium so that the total force 
on each particle vanishes = 0. Then the scalar product 

Ft-Sri (2.14) 

which is the virtual work of the force F in the displacement 6f also 
vanishes. Hence the sum over all bodies 


££ -6fi = 0 (2.15) 


which so far says nothing new. 

If we now decompose the force into its external component F ext , 
and internal (constraint) component F int , 


Fi 


F‘ xi 



(2.16) 


we get that 

• 6 ? i) = 0 ( 2 - 17 ) 

t 

Suppose we now restrict ourselves to systems where the the net virtual 
work of the forces of constraint is zero 7 . Then our equation 2.16 
becomes 

£F‘ xt .6r i =0 (2.18) 

« 

7 This is true for rigid bodies and will be assumed for all systems considered from 
now on. It can fail if, for example, rolling friction is present. 
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f zx 


Figure 2.5: Virtual work principle - the lever. 


The equation 2.18 is the static principle of virtual work. 


Example of static principle of virtual work: Consider a lever with one 
arm length a and the other b as shown in Figure 2.5. External forces A and 
B are applied at each end. Let x and y be the heights of each end. Since the 
lever is rigid and so has only one degree-of-freedom, we have a constraint 
between x and y such that x/a = y/b. The principle of virtual work states 
that Adz + BSy = 0 for kinematically admissible displacements, i.e. ones 
for which 8y = -b8x/a . Thus the principle gives that (A — bB/a)8x = 0 
for every virtual displacement 6x and so the condition for equilibrium of the 
lever is Aa = Bb y i.e. the moments must balance. Note that if we did this 
problem using Newton’s second law, we would have to explicitly include the 
reaction force from the fulcrum, in addition to the external forces A and B. 
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2.2.3.4 D’Alembert’s Principle 

In order to obtain a principle which is valid for more general motion, 
we use a device that originated with J. Bernoulli and developed by 
D’Alembert. Start with the equation of motion equivalent to 2.1: 

Fi = ft (2.19) 

This can be written 

Fi-pi = 0 ( 2 . 20 ) 

which can be interpreted as saying that the particle is in equilibrium 
with the total effective force Fi — p, made up of external and inertial 
parts. Thus instead of 2.18, we have 8 

£( ft"' - ft) Sr, = 0 (2.21) 

t 

This is called D’Alembert’s principle. 


Example of dynamic virtual work principle: Consider a pendulum of 
mass m supported by a chord of length r. See Figure 2.6. It has gravity 
as its only external force. The virtual work equation is 8x r (mx+ tngz) = 
0. But x = (rsin(ff),r(l — cos(0))) where 9 is the angle the chord makes 
with the vertical. By the chain-rule, 8x = (rcos(0),-rsin(0))50 and our 
virtual work equations becomes 89(9-\-g sin(0)) = 0 Since the virtual rotation 
89 is arbitrary, the dynamic equation of motion for the pendulum is 9 + 
(<7/r)sin(0) = 0. This is the classical equation for the pendulum. Notice 
that any dependence on the mass has disappeared. It can be solved exactly 
in terms of elliptic functions, or approximated by the harmonic oscillator 
equation for small 9 . 


2.2.3.5 Kinematically Admissible Virtual Displacements 

A virtual displacement is kinematicaHy admissible when it is con¬ 
sistent with the constraint conditions $ = 0. By this we mean that it 


8 Again, we are assuming that the forces of constraint produce no virtual work. 
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Figure 2.6: Virtual work principle - the pendulum. 

satisfies the variational form of the constraint condition 

$ f Sr = 0 (2.22) 

The equation 2.22 follows directly from the chain rule since 6 acts like 
differentiation 9 . 

We can interpret equation 2.22 geometrically as follows: The con¬ 
straint equation defines a hypersurface in phase space 10 . At each 
point of the surface, there is a tangent hyperplane. The rows of the 
constraint Jacobian form a basis which spans a complementary sub¬ 
space to the tangent space. For example, if we have a space with N 
coordinates and the row rank of the constraint Jacobian is rk, then 
the tangent space is of dimension N — rk. Thus the row rank of 
determines the the dimensionality of the hyperplane, and hence the 
constraint surface. Equation 2.22 simply says that the kinematically 
admissible displacements are those that are orthogonal to this com- 

9 This can be shown by taking Taylor’s expansion of 4*(r + 6r) at r to first order, 
way of stating that the solution to an equation will define some higher di¬ 
mensional surface in the coordinate space. 
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plementary space. Thus the admissible displacements are those living 
in the tangent hyperplane. 


A good example is when we have a single constraint equation $(r) = 0. 
In this case, the constraint Jacobian is just the transpose of the constraint’s 
gradient V$. But this gradient is always perpendicular to the surface 11 and 
the vectors which are perpendicular to it are by definition the elements of 
the surfaces tangent plane at that point. So the kinematically admissible 
displacements are the elements of the surfaces tangent hyperplane. 


2.2.3.6 An Essential Difficulty 

We must now make a crucial observation: We would like to take Equa¬ 
tion 2.21 and claim that each term in the sum must vanish individ¬ 
ually. This cannot be done as the equation stands because Equation 
2.21 is in terms of the underlying Cartesian virtual displacement’s 
6fi which must still satisfy all admissible kinematic constraint condi¬ 
tions the system is subjected to 12 . We cannot conclude because the 
total vector sum must vanish, or even that each individual term must 
vanish, that for each » 


{Ft*-fr)-6* = 0 (2.23) 

We could do this if the individual 6fi were without constraint, since 
then they are arbitrary, not having to satisfy and admissibility con¬ 
ditions. Imposition of constraints on the system’s motion produces 
relations between the various components of each 6r, which destroy 
the required independence. A very fancy way to say this is that the 
principle of virtual work requires only that the vectors 

Ft” ~ P, 

11 This is the general property of gradients. V/ is always perpendicular to any 
level surface f = const . 

12 In geometric terms, they are only elements of the lower dimensional tangent 
spaces to the constraint surface. 
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be orthogonal to the subspace of virtual displacements that are ad¬ 
missible, and so, in geometric terms, are elements of the constraint 
surface tangent space. 

The way around this problem is to introduce either generalized coor¬ 
dinates which have the constraints built in and so the virtual change 
in the coordinates are independent, or to introduce what are called 
Lagrange multipliers which also allow for independent virtual dis¬ 
placements, but at the cost of adding an extra unknown vector. We 
now review generalized coordinates and then discuss the variational 
equations of motion. 

2.2.4 Generalized Coordinates 

Sometimes it is easy to pick a new set of coordinates that are more 
natural to the problem at hand 13 and have all or most of the system 
constraints already built into them. These are then the generalized 
coordinates 14 for the system. 


Examples of generalized coordinate systems: 

1. Pendulum - angle between attachment chord and vertical. 

2. Planetary motion about the sun - polar coordinates in plane of orbit. 

3. Macpherson front suspension - vertical translational “jounce” and 
steer angle of spindle. 


What we require for generalized coordinates is that we can write 
the Cartesian coordinates for each body uniquely in terms of the 
generalized coordinates 9 x,<? 2 ) • • • ,<7 n . Thus the time variation of the 
generalized coordinates completely determines the dynamics. The 
minimal number of generalized coordinates needed in describing the 
system is called the degrees of freedom 15 of the system. We will 

13 By this we mean they may considerably simplify the problems formulation. 

14 Which need not have the dimensions of length 

16 Written often as DOF’s. 
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see later how to know when one has a correct set of generalized co¬ 
ordinates given a set of constraint equations, but it should be clear 
to the reader that generalized coordinates must number at least N 
- rk(€f) K , so that enough coordinates are available to describe the 
motion. If more are given, then they are not independent, and must 
satisfy constraint conditions. 

Examples of counting DOF’s: 

• One DOF: 

1. A piston moving up and down. A translational joint. 

2. A rigid body rotating about a fixed axis. A revolute joint. 

• Two DOF: 

1. A particle moving on a surface. 

• Three DOF: 

1. A particle moving in space. 

2. A rigid body rotating about a fixed point (a top). 

3. A spherical joint. 

• Four DOF: 

1. Two components of a double star revolving in the same plane. 

• Five DOF: 

1. Two particles at a constant distance form each other. 

• Six DOF: 

1. Two planets revolving about the sun. 

2. A rigid body moving freely in space. 


Example of constrained particle cases: 

10 Where rk(M) is the rank of the matrix M. A theorem in linear algebra states 
that row rank = column rank = rank of a matrix. 
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1. If a particle were constrained to move on a curve, we would need only 
a single generalized coordinate. 

2. If, for example, the particles were constrained to move on a plane, 
there are fewer independent coordinates needed to completely describe 
the system. Each particle now requires only two coordinates, so the 
total system requires 2M. We have imposed M independent constraints 
on the particles of the system. Thus only 2M generalized coordinates 
would be required to fully describe the motion. 

3. Suppose in the general case we have M particles and have imposed k 
constraints. Then as in section 2.2.2 we would intuitively think that 
we need only q\, 92 >... , 93 M-k generalized coordinates if the Jacobian 
has maximum row rank, where we have 

r 2 = rl(qi,q 2 , ..,qsM-k,t) ^ ^ 

= fl(q it q2,- ■ - ,qsM-k,t) 

The number SM-k is the total number of degrees of freedom for the 
system. 


2.2.4.1 Change of Variables to Generalized Coordinates 


Suppose we have our M bodies described by a set of N generalized 
coordinates. We always assume that the relation Equation 2.24 is 
differentiable, so that we can relate the derivatives of the two systems 
of coordinates. This is done using the chain rule and the Jacobian 
of the coordinate transformation. The Jacobian is the matrix fq 
of partial derivatives 


gfl .1 

dr i.j 

dqi 

is* 

3r,.j 


/•*> 
dr i.. 


dqi 

dqn 


(2.25) 


. . . dfMj* I 

dqi dqs ) 


\ 
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Example of a single particle Jacobian: If we have only one particle whose 
coordinate vector is r, then its Jacobian matrix with respect to the m gen¬ 
eralized coordinates 7i, 92 , • • •, q m will be 


r 4.‘.J 



(2.26) 


where 1 < : < 3 and 1 < j < m. 


2.2.4.2 Chain Rule 

We need the general version of the chain rule. We just state the result 
and leave the proof to the references. The chain rule tells you how 
to change variables when you are dilferentiating. Suppose y = y(x), 
then 17 


3. \ 

f 

9 _l 


dvi 

— 

a 


V dy n / 

V 


or, in much neater matrix notation, 


d»A 

... 

Ml 

dZ\ 

8x2 

dx m 


dy* . .. 

ft. 

dx\ 

6x2 

dXm 

Sj/n 

dyn . . . 

#.Vn. 

dz 1 

dX2 

dZm 

notation, 



( 3- 
a V_ 

dX2 


\ 


J \ aZZ J 


(2.27) 


A = „_A 

dy Vz dx 


(2.28) 


As a special case, we have the time derivative of a vector with 
respect to the time derivatives of the new coordinates 


y = VxX 


(2.29) 


Example of a particle on a curve: A position on a space curve is deter¬ 
mined by a single generalized coordinate q 18 . Then the particles Cartesian 
position r can be written as 

r = r{q) (2.30) 

A7 We use operator notation in the next equation, since the operand is a dummy 
place-holder anyway. 

18 Which could be taken as arclength along the curve. 
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After differentiating, the particle’s velocity is, using the chain rule, just 

r = r q q (2.31) 

where the Jacobian r q is the curves tangent vector. The particles acceler¬ 
ation vector is then 

r = r q q + r qq q 2 (2.32) 

The first term on the right can be interpreted as the generalized acceleration 
q and the second term as a centripetal acceleration q 2 . This is obvious in 
the case of uniform circular motion, since then q = 0. 


Example of a particle on a surface: In this case we have two generalized 


coordinates q = [gi,fl2) r > 

and the coordinates are r = f(q). 

The velocity 

relation is now 


r- r q q 

(2.33) 

and acceleration 


r = r q q+ {r q q)$q 

(2.34) 


Again, the first term is the generalized acceleration, while the second term, 
which is the Jacobian of a vector derived form the Jacobian, is a combination 
of mysterious accelerations such as the coriolis and centripetal. The form of 
this equation remains the same for any general set of coordinates. 


Example of an application to holonomic constraints: If we take a holo- 
nomic constraint and differentiate, we get 

(2.35) 

which gives a condition that the velocities r must satisfy in order to be 
consistent with the constraint. If the constraint is time independent, then 
the right hand side of the equation will vanish, and so the constraint Jacobian 
will annihilate the velocity vector, that is to say, the vector is in the null 
space or kernel of the Jacobian matrix. 

If we differentiate the constraint one more time, then we have 

$ f ?=-7 (2-36) 


where 


7 = ('M')^ + 2 ®rt + 


(2.37) 
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If we work instead with generalized coordinates, we have 

$<f = $ frq 

and the equations 2.35 and 2.36 become 

$tf = 

$tf = -7 

where now 

7 = (M)j? + 2 $qt + $tt 


(2.38) 


(2.39) 

(2.40) 


2.2.4.3 Kinematically Admissible Generalized Virtual Dis¬ 
placements 

In the same way as section 2.2.3.5, we have kinematically admissible 
generalized virtual displacements 6q. Again, by the chain rule we 
have 

$tfq = 0 (2.41) 

where now 

2.2.4.4 D’Alembert’s Principle in Generalized Coordinates 

We now transform the virtual work Equation 2.21 to generalized co¬ 
ordinates and forces. 


First, we consider a single body with generalized coordinates q. From 
r = r(g), we have the generalized virtual displacement 


c* 

-u 

II 

(2-42) 

and the time derivatives 


r = rtf 

(2.43) 

and acceleration 


r = rtf + {uq)tf 

(2.44) 
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We then can write equation 2.21 in terms of all kinematically admis¬ 
sible 6q as 

+ (?}?),?) - SW* = 0 (2.45) 

where 6W ext is the external virtual work given by 

6W ext = 6fr?F ezt (2.46) 

The quantity Q ext = r?F ezl is called the generalized external force. 
After collecting terms, we can write equation 2.45 as 

6q r (mr~r^q + mrj- Q ext ) = 0 ( 2 - 47 ) 

for all kinematically admissible virtual generalized displacements. 


Thus if we have an independent set of generalized coordinates, all 
virtual displacements are kinematically admissible and we must have 

mrT r^q + mf? {r^q)^q - Q ext = 0 (2.48) 

Note that by picking independent generalized coordinates, we have 
solved our essential difficulty! This is because we are no longer re¬ 
stricted to some subspace of virtual displacements which must be 
orthogonal to the vector 

(2.49) 

but now every vector in the space of virtual displacements must be 
orthogonal to 2.49 and so the vector must vanish. 


The case of multiple bodies is identical. One just must add an index 
* to each of the quantities above. The crucial observation is that one 
may choose an independent set of generalized coordinates. This gives 
us 


£ miffifttf - Q ext + £ mSi [ (f if q)jq = 0 


(2.50) 


* 

where the summation is over 
generalized force is 

Q ext 


i 

the number of bodies, 




-* ext 


and the external 


(2.51) 
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2.2.5 Lagrange Multipliers 

In this section, we show an alternative way of dealing with deriving 
the equations of motion with respect to constrained coordinates. With 
generalized coordinates, we eliminated the constraints by choosing a 
proper Bet of new independent coordinates. With Lagrange multipli¬ 
ers, the constraints are dealt with directly, and we incorporated into 
the equations of motion. 


2 .2.5.1 Lagrange Multiplier Theorem 


Haug in reference [37] calls the following theorem the Lagrange mul¬ 
tiplier theorem. Suppose we have b in R n and A a m x n matrix. 


If 


5^6 = 0 


(2.52) 


for all s in R n such that 


As = 0 


(2.53) 


then there exists a vector A in R m , called a Lagrange multiplier, 
such that 

s r b + f'A T X = 0 (2.54) 

for all s in R n , and if A has full row rank, then A is unique. 

It is worth looking at this from a geometric point of view. We 
claim that it says something fairly interesting geometrically. The first 
condition, Equation 2.52, can be interpreted as saying 6 is orthogo¬ 
nal to the subspace of vectors that A takes to zero (equation 2.53), 
b±ker(A) 19 , and the conclusion, equation 2.54, is that b is the image 
of some vector —A under the matrix transformation A T . This last 
statement can be written 6 eim(A T ). That either this is true or the 
first condition must then be false, and there is a vector with s^b 7 ^ 0 , 
is called the Fredholm alternative theorem. It follows essentially 
from the fact that if a vector is in the kernel of a transformation, there 
is no way to recover it as the inverse image under the transformation 
of some other vector. 

10 As we have said before, the set of vectors s such that As = 0 is called the kerntl 
of A. 
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2.2.5.2 Application of the Theorem 

Again, we will restrict ourselves to a single particle with coordinates 
r for simplicity. 

We know that 

6r r (m?-F)=0 (2.55) 

for all admissible virtual displacements 6r that satisfy 

= 0 (2.56) 

Now, in the Multiplier theorem, take 

b = mr — F 


and 

A = $ r - 

Then we must have a A such that 

6r r (mF - F + = 0 (2.57) 

for all 6f. We have done what we wanted to, we now have an equation 
that must hold for every virtual displacement, and so the coefficient 
of 6r must vanish. Our Lagrange Multiplier form of the con¬ 
strained Cartesian dynamics equations are then 

(mr — F + A) = 0 (2.58) 


It will be very important for numerical applications to note that equa¬ 
tion 2.58 is not an ordinary differential equation, but, because of A, 
is instead an differential-algebraic equation. These demand special 
numerical techniques for their solution. 

For multiple particles, we write 


-T }T 


' _ =fT zfT -fT 

~ l r l > r 2 i ' * * > r N) 


r _ I jpcztT jptztT j?eztT\T 

* “ 1*1 1 * ## 











\ 






and 
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for the particle position vectors and external forces. The mass matrix 
is then the diagonal matrix M = diag(mil, • • •, niff I) 20 . The multiple 
particle version of equation 2.58 is then 

[M~r-F + $T\)=0 (2.59) 

In the generalized coordinate formulation using Lagrange multi¬ 
pliers 21 , we use the Jacobian as our coordinate transformation. If we 
proceed as above, but use the chain rule for taking derivatives with 
respect to our generalized coordinates, we get the dynamical equations 

+ rjAf (ftf)tf - Q = 0 (2.60) 

2. 2. 5.3 Matrix Formulation of Constrained Equations of 
Motion 

We can combine all our equations into one nice matrix equation. By 
combining equations 2.36 and 2.59, we get 

(s:7)GM7) 

Haug in reference [37] shows that this equation, when combined with 
the obvious proper initial conditions, provides a complete description 
of the system dynamics as long as M is positive definite for all kine¬ 
matically admissible vectors 22 

Similarly, we can form the matrix version of the generalized coor¬ 
dinate constrained dynamics equations using equations 2.39 and 2.60. 
First, it’s useful to rewrite equation 2.60 as 

M'q + $T\ = Q“‘ (2.62) 

where we have used the reduced mass matrix 

M = f[Mf f (2.63) 

20 Where I is the identity matrix of the proper dimension. 

21 We can indeed choose generalized coordinates that themselves are subject to 
additional constraints. 

22 A matrix M is positive definite if for all vectors a, a T Ma > 0. This is 
equivalent to having only positive eigenvalues. 
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We also define a new force that combines our generalized force with the 
pseudo-force terms that arise from the non-inertial reference frames 
of the generalized coordinates 

R = Q‘" - (2.64) 


We then have as our matrix equation 



(2.65) 


Haug [37] shows that this is a unique prescription for the system’s 
dynamics whenever the generalized coordinates are non-singular and 
the constraint Jacobian has full row-rank. 
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Figure 2.7: Rigid Body. 

2.3 Rigid Body Kinematics 

This section reviews rigid body kinematics. It introduces both Euler 
angles and Euier parameters. In particular, it gives the description of 
rigid body motion required by the vehicle dynamics model. 


2.3.1 The Rigid Body 

We describe some of the essential aspects of rigid bodies. 

2.3.1.1 Definition of a Rigid Body 

We think of “solid” bodies as objects rigid in their shape, with each 
point of the body holding a fixed position with respect to all the other 
points. More precisely, a rigid body is a large collection (thought of 
as infinite in number) of particles labeled by their position vectors r,, 
and which have a constant distance aj = u — r ; ) 2 between each pair 
of particles. See Figure 2.7. As the rigid body undergoes acceleration, 
reaction forces must exist in order to hold the rigid body together, 
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satisfying the holonomic constraint conditions. The constraint con¬ 
ditions of fixed distances between all points insures that all angles 
between particles are also fixed. 

A rigid body is the first approximation to most engineering struc¬ 
tures, when we can ignore the compliance of the system, and treat it a 
a perfectly rigid object. Such an object has no internal energy 23 since 
the internal constraints do no work. If internal energy is important 
for the problem at hand, such as when the energy stored in a spring 
is a significant part of the total system energy, then the rigid approx¬ 
imation fails. The present dynamics models treat many parts of the 
vehicle as rigid bodies, since the elastic energy for those parts is of no 
importance to the total system. Thus the strain energy of, say, the 
front Swiss arm is probably not crucial to the dynamics, but the coil 
springs must be treated as compliant bodies. 

We remember that a rigid body is fully described by sti indepen¬ 
dent coordinates or degrees of freedom; we must specify the position 
of a single point within the body (usually its center of mass. Cf. the 
next section) and its orientation. Both the point and the orientation 
require three parameters to be determined in space. We’ll quickly 
review the center of mass, and then take some time thinking about 
the idea of how a body’s orientation is best described. 

2.3.1.2 Center of Mass 

The center of mass is the position R at which all the inertial mo¬ 
ments vanish, and would be the point (not necessarily in the body) at 
which if it where supported there, it would not rotate if the support 
were accelerated. We can calculate the center of mass R for a set of 
particles as 


and, for a solid rigid body, as 

(2.66) 

S fm(x)r{x) 
fm(x) 

(2.67) 


where the integrals are taken over the solid’s volume. 


23 It has, for example, no elastic or plastic strain energy. 
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The center of mass turns out to be a particularly good place to take 
as the origin of coordinates in many dynamics problems. 


Examples of centers of masses: 


1. The center-of-mass of a disk with uniform density is at its geometric 
center. 


2. A front engined car chassis with engine has a center-of-mass forward 
of its geometric center due to the mass of the engine. The exact c.g. 
shifts under vehicle loading. 


2.3.1.3 Rigid Body Orientation 

The angular orientation of the rigid body is measured by how an 
orthogonal coordinate system fixed in the body 24 changes its angular 
relationship to a orthogonal coordinate system fixed in the global iner¬ 
tial system 26 (see 2.8). We can measure this by the direction cosines 
of the body fixed coordinates with respect to the world coordinates. 
If we call the unit vectors in the £*- system 



and in the x-system 

i,j,k 


24 We will denote this as the x'-system. 

25 Denoted as the x-system. 
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ai = a •» 
a 2 = £ • J 

0-3 = i 1 • k 

0i = 

02 = f‘l ( 2 - 68 ) 

03 = f ■ ^ 

71 = & •» 

72 = *' • 7 

73 = ft ■ k 

These are called direction cosmes since the dot product of unit vectors 
is just the cosine of the angle between them. The direction cosines 
are the components of the x* unit vectors with respect to the x unit 
vectors. 

Because each of the set of coordinate unit vectors are orthogonal , that 
is, they are pairwise at right angles with one another for our Cartesian 
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coordinate systems and are of unit length 26 , they form the components 
of an orthogonal matrix A 27 

( «i fix 7i \ 

«2 A 72 (2.69) 

<*3 /?3 73 ) 

We could describe the orientation of the body by the direction cosines, 

but we would need nine parameters connected by the six conditions 
of forming an orthogonal matrix. This would be very inefficient, but 
could be done. Instead, we will try to find a smaller set of parameters 
that will also do the job, but first we review rotations in Cartesian 
three-space. 

2.3.2 The Description of Rotations 

We now describe the theory of rotations in three-space. This can be 
made quite abstract, but here I only present a utility grade theory. 
The most important fact about rotations is that you can’t add them 
in the same way as you would vectors. They are different beasts 28 , 
and only infinitesimal rotations are true vectors in space. Only the 

2G These conditions are of the form i* ■ j f — 0 and k* • k* = 1, etc. 

27 An orthogonal matrix is one whose transpose is its inverse, so A T = A~ l . This 
is equivalent to the columns being orthogonal unit vectors. 

28 They are in fact elements of a non-commutative Lie group of real 3x3 orthogonal 
matrices, with unit determinants, called SO(3). They naturally form a smooth 
manifold which can be locally parameterized with three coordinates in a manner 
such that all group operations are also smooth functions of the parameters. The 
Lie group SO(3) is not simply connected, but has another group called SU(2) of 
complex 2x2 unitary matrices with unit determinant which is its simply connected 
covering space. The Lie algebra associated with SU(2) are the 2 Hermitian matrices. 
The real form of these, the real 3x3 skew-symmetric matrices form the Lie algebra 
of SO(3). The infinitesimal rotations are the members of this Lie algebra of SO(3), 
and so can be associated with 3-dimensional vectors. In real terms, SU(2) is the unit 
sphere in real 4-space. The standard real parameterization of this sphere will be the 
four Euler parameters e 0l e x » c 2i c 3 W ^ 1 normalization e 2 + e 2 4* t\ 4* e 2 = 1 that 
I will introduce in a much more elementary manner in Section 2.3.2.4 (a complex 
parameterization by two Caley-Klein parameters is also sometimes used). The Lie 
algebras of each map to their corresponding group via the matrix exponential. A 
basis for the Lie algebra of SU(2) is the identity matrix and the three Pauli spin 
matrices <r t *. The two Lie algebras can be identified by the mapping x ■—♦ x • o, 
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definition of Euler angles and Euler parameters, along with the finite 
and infinitesimal rotation formulas are critical. 


2.3.2.1 Rotations as Orthogonal Matrices 

The reader should be familiar with the fact that the orthogonal matrix 


M t {B) = 


f cos(0) 
— sin(0) 
0 


sin(0) 0 \ 
cos(0) 0 I 
0 l) 


(2.70) 


is a counter-clockwise (right-handed) rotation about the z-axis of 9 
radians. Similarly, we can write down rotations about the y-axis 


and x-axis 



( cos(0) 0 

sin(0) 

m,w = 

0 

1 

0 


- sin(0) 0 

cos(0) 


(i 

0 

0 

M x (9) = 

0 cos(0) 

sin(0) 


^ 0 - sin(0) 

cos(0) 


(2.71) 


(2.72) 


A crucial (and complicating property) of rotations is that the order 
in which you perform finite rotations matters 29 ! The reader is urged 
to check that a rotation of a body about the x-axis followed by a 
rotation about the y-axis is not in general the same total rotation as 
if we first rotated about the y-axis and then rotated about the x-axis. 
M x M y ^ MyM x , so M z and M v do not commute. 


2.3.2.2 Euler Angles 

As noted, the nine direction cosines are not a good set of orientation 
coordinates since they are subject to six relations. To find the min¬ 
imal independent set (by dimension counting, there are three since 
the six relations are independent) of three parameters, we break the 

while the map between SU(2) and SO(3) is a 2-fold covering transformation and so 
half-angles will start to appear in the formulas. 

29 We will see later that for an infinitesimal rotation order does not matter since 
they are true vectors, and so are commutative in their addition. 





Figure 2.9: Diagram illustrating the three successive rotations 
through Euler angles a, /?, and 7 . 


angular motion into three rotations about special axes. These will be 
the Euler angles. 

An arbitrary orientation of the coordinate axis (and so of the rigid 
body) can be obtained by the following sequence of successive rota¬ 
tions about single axes 30 

1 . Rotating (xyz) through an angle a about the z-axis to give 
( x'y'z '). The orthogonal matrix is given by M z (a). 

2 . Rotating (x'y'z 1 ) through an angle /3 about the y'-axis to give 
(x"y"z"). The orthogonal matrix is given by M y < (/?). 

3. Rotating (x"y"z") through an angle 7 about the 2 "-axis to give 
( x'"y'"z '"). The orthogonal matrix is given by M z »( 7 ). 


This is illustrated in Figure 2.9. 

30 Remember, the order in which we perform the rotations matters! 
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(2.73) 


The complete orthogonal transformation is 

M(a,(3, 7 ) = M z {a)M v >{0)M t n( 7 ) 
and so M(a,/?, / y) = 

1 cos(/J) cos(a) cos(- 7 ) cos(/3) sin(a) cos( 7 ) - sin(/?) cos( 7 ) ^ 

- stn(a) sin('y) cos(a) sin( 7 ) 

-cos(/ 9 )cos(a)sin( 7 ) -cos(/?)sin(a)sin( 7 ) sin( 0 )sin( 7 ) 
-sin (a) cos ( 7 ) cos(a)cos( 7 ) 


sin(/?) cos (a) 


sin (/?) sin(a) 


cos(/3) 


(2.74) 

We must warn the reader that there is nothing special about which 
axes the rotations are to take place. There are other conventions in 
use. But the sequence given above is the standard modern physics 
one. An alternative that occurs in engineering is to use a yaw, pitch , 
and roll sequence. 

The Euler angles are a minimal set of rotation coordinates. They 
can be used to find orientation from Euler angle rates since they are 
additive. But the reader is encouraged to check that the transfor¬ 
mation Jacobian of Equation 2.74 is singular at certain values of the 
parameters, thus it is not a good set of coordinates for numerical cal¬ 
culations, since extra care would be involved in avoiding the singular 
nature of the coordinates. A better set of coordinates are the Euler 
parameters. But first, we motivate their existence by deriving what is 
called the finite rotation formula. 


2.3.2.3 Finite Rotation Formula 

Suppose we are given a unit vector rf, another vector r, and ask for the 
new vector r* which arises when we rotate r through an angle $ about 
the axis vector rf. The answer is elementary and can be obtained from 
geometry by considering Figure 2 . 10 . Then 

r* = ON + NV + VQ 

= ijin-r) + (r- rf(f} -r))cos($) + (if x r)sin($) (2.75) 





Figure 2.10: Derivation of iinite rotation formula. 


This gives 

r* = r cos($) + tj(ii • r)(l — cos($)) + (rj x r) sin($) (2.76) 
This formula works for any magnitude of rotation. 

2.3.2.4 Euler Parameters 

If we start from the finite rotation formula Equation 2.76, there is 
a magic parameterization of it in terms of four Euler parameters 
e 0 ,c = {ei,€ 2 ,e 3 ). These are related to a finite rotation by angle x 
about an axis unit vector f( by 


e 0 = cos(x/2) 

e = tjsin(x/2) (2.77) 

Euler parameters clearly satisfying a normalization condition e* + 

e i "h e 2 "h e i = 

Then from trigonometric identities, 

cos($) = el - e\ - t\ - t\ 


(2.78) 
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and 

rf sin($) = 2coc (2.79) 

Equation 2.76 and then be written 

7 = r(e* - e\ - e\ - e’) + 2e(e • r + 2(e x r)e 0 (2.80) 

If we now expand Equation 2.80 component-by-component, we can 
recover the cosine matrix M which relates r* to r. It is 

( e o + e i ~ e l - e l 2 ( e i c 2 - e o e s) 2 (e 1 e 3 + e 0 e 2 ) 

M = I 2(eie 2 + e 0 e 3 ) 4 - e\ + e\ - e\ 2(e 3 e 3 - e 0 ei) 

* ‘ 4 


^ 2(eie 3 + € 0 ^ 2 ) 2(e 2 e 3 — cqC 1 ) 4 c\ t\ + e 


(2.81) 


If we define 



( -ei 

Co 

—e 3 

e 2 

E = 

C2 

C3 

co 

-Cl 


-e 3 

-e 2 

Cl 

e 0 


( —e\ 

eo 

C3 

-e 2 

G = 

-e 2 

— e 3 

e 0 

Cl 


" e 3 

e 2 

“Cl 

e 0 


(2.82) 


(2.83) 


then 


M = EG t (2.84) 

We can then relate the Euler angles to the Euler parameters by 
the formulas 


,4> + 0 ^ ft \ 

e 0 = cos(^—)cos(-j 

<h — tb . ,9 v 

Cl = cos(^y^-)sin(-) 

e 2 = sm(—)sin( 2 ) 

. ,4> + \ 

e 3 = sin(—)cos( 2 ) 


(2.85) 


The explicit form of this construction is very important, since it 
provides a way of parameterizing rotations in a non-singular manner. 
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The reader is encouraged to check that the Jacobian of the mapping 
is everywhere of full rank as promised. It is in terms of the Euler 
parameters that the numerical solution of rigid body motion can be 
performed most readily since the parameters are always good and so 
provide a stable method of describing the motion. 


2.3.2.5 Infinitesimal Rotations 

Suppose if instead of a finite angle $, the rotation is infinitesimal in 
the sense that it is so small we can always ignore second- or higher- 
order terms in any Taylor series expansion involving $. Thus is we 
make the approximations cos($) « 1 and sin($) « $, Equation 2.76 
becomes 

r* = r + ${r}X r) (2.86) 

We can use Equation 2.86 to write the effect of a infinitesimal 
rotation on a vector r in the concise form 


Af($)r = Ir + x r 
= (I + $7?x)r 

= (/-$x)r (2.87) 


where $ = $rj. 

A useful device is to notice that the operation of “taking the cross 
product” ox is the same as multiplication by the skew-symmetric 
matrix a generated by the vector a 


( 0 — a 3 


a = 


o 3 
V ~ a 2 


0 

a t 


0-2 ^ 
—ai 

0 ; 


(2.88) 


It is easy to check that axx = ax. In particular, we can replace x 
by $fj. 


2.3.3 Rate of Change of a Vector 

In this section, we find the kinematic description of velocity and ac¬ 
celeration for a rigid body rotating with respect to a fixed inertial 
frame. Rotation angles will be denoted 9, and the angular velocity Q. 
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2.3.3.1 Infinitesimal Rotational Relations 

From equation 2.86, we see that for a small rotation 0, we can write 
the rotation’s affect on a vector r as 

M(0)r = r + Or (2.89) 

Thus the infinitesimal change dr is 

dr = M{6)r-r 

= Or (2.90) 

Because we want to require that 0 be infinitesimal, we’ll write it dO = 

cD 31 ,so 

dr = dOr (2.91) 


2.3.3.2 Vector Rate of Change Equation 

We see that the rate of change, or velocity, due to a rotation alone, is 

dr/dt = <Dr (2.92) 

It should be clear to the reader that the total rate of change of any 
vector observed in a fixed inertial reference frame , is the total of the 
change as observed in a frame fixed in the body and the change due 
to the relative rotation rate w of the two frames. Thus the total rate 
of change of any vector G 

[dG / di)gpa Ct = ( dG/dt)t 0 dv + <DG (2.93) 

2.3.3.3 Velocity and Acceleration 

If we now take G to be the position vector of a point on the rigid 
body, we have 

(dr/dt) !pace = (dr/dt) b 0< i v + (Dr (2.94) 

and so in terms of the velocities 

Vspace = Vbody "I" (2.95) 


31 Remember, the infinitesimal rotations, in contrast to finite rotations, d really 
are vectors, and so can be added. 
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If we apply our equation 2.93 to v gpace , we get a famous and important 
result relating the accelerations in the two frames 

{dv space / dt) tpace = ®«pac« 

= {dV)p aee jdt)i,odfi space (2.96) 

= a body + 2 (u>v bodv ) + wwf 

Equations 2.95 and 2.96 are our fundamental kinematic equations for 
rigid bodies. 

The two “extra” terms in the right-hand side of equation 2.96 are 
of special interest and should be known to the reader. They exist 
because of the non-inertial rotating frame. The term 

dcorioli, = 2(tiv iodv ) (2.97) 

is called the Coriolis acceleration, and h 

a,centrifugal = WWr (2.98) 

is called the Centrifugal acceleration. The Coriolis force exists only 
when motion is occurring with respect to the rotating frame. The 
angular velocity enters as the first power, while the centrifugal force 
depends upon the angular velocity to the second power. 

2.3.4 Euler’s and Chasles’ Theorems 

Two theorems about rigid bodies are relatively easy to establish by 
considering the eigenvalues of an orthogonal matrix. Any 3x3 or¬ 
thogonal matrix has only one real eigenvalue A = 1, the other two are 
complex conjugates with non-zero imaginary part. The eigenvector 
corresponding to A = 1 is an axis vector for the rotation, while the 
other’s are complex. 

Euler’s theorem on the motion of a rigid body states that 

Theorem 1 (Euler) The general displacement of o rigid body with 
one point fixed is a rotation about some axis. 
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Because any vector along the axis of rotation is left fixed, it must 
be an eigenvector of the transformation, with a unit eigenvalue. The 
proof of Euler’s theorem consists of showing such an eigenvalue exists 
for any transformation with a fixed point. 

A corollary of Euler’s theorem is called Chasles’ theorem which 
states that 

Theorem 2 (Chasles) The most general displacement of a rigid body 
is a translation plus a rotation. 

These are the final things we have to say about rigid body kine¬ 
matics. 
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2.4 Dynamics of Rigid Bodies 

In this section, we review the inertia matrix and present the Newton- 
Euler equations of motion for a rigid body. A theorem already men¬ 
tioned due to Chasles’ states that an arbitrary motion of a rigid body 
can be represented by a translation plus a rotation. This suggests that 
we can decompose the dynamics of a rigid body into two independent 
parts: Translation and rotation. The translation part is easy and just 
involves the center-of-mass motion. The rotational part involves the 
concept of angular momentum and the inertia matrix or tensor. 

2.4.1 Rigid Body Motion 

We can treat a rigid body as a collection of particles held together by 
some unknown (and arbitrary) set of internal forces which act along 
a a line between the particles in an “opposite and equal” fashion. We 
know that the translational and angular momentum are conserved in 
the absence of external forces and torques, and if, indeed, such forces 
and torques are present, we have the dynamical equations 

F ^ext ^ 

- p 

and 

N exl = L (2.100) 

We have also previously remarked that the angular momentum de¬ 
pends upon the choice of origin of coordinates. 

If we have two coordinate systems attached to the rigid body as 
shown in Figure 2.11, then an essential property of the body is that 
all points move together, translating and rotating at the same rate. 
Thus we have both angular velocity vectors for the two coordinate 
systems being equal and arrive at the conclusion that the angular 
velocity vector of a rigid body, with respect to any coordinate system, 
are equal 32 . We can prove the following relationships 

t>2 = V\ 4* CjR 

U>2 = Wi 


(2.99) 


32 The proof of this is easy, and can be found in |30). 


( 2 . 101 ) 

( 2 . 102 ) 
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Figure 2.11: Relation between sets of rigid body coordinates with 
different origins. 

Note that translational motion does not have this absolute relation¬ 
ship. This relation implies, when combined with Equation 2.95, that 
at any instant, a new center of coordinates can always be chosen so, 
with respect to this new coordinate system, the body undergoes a 
pure rotational motion, without any translation. 

2.4.2 Angular Momentum 

The total angular momentum L of a rigid body is 

L = j v {rxp)dV (2.103) 

Its value depends upon from where r is measured, that is, form the 
choice of the origin. If we pick a second origin, so that r = r' + a, 
then 

L = L' + axp (2.104) 

We see that the angular momentum depends on the choice of origin, 
except when the system is a rest as a whole, so that p = 0. If we now 
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take as our coordinate system the center-of-mass frame, then we can 
replace a = R. 

2.4.3 Inertia Tensor 

If we consider how to express the kinetic energy of a rigid body, we 
find that there should be a linear relation between the body’s angular 
momentum and its rotational velocity. We can write it as 

L = JQ (2.105) 

where J is the inertia matrix (or tensor since it is a doubly-indexed 
quantity J it j with the proper transformation property for a second 
rank covariant tensor). 

We can write the inertia tensor as 

Jij = J v p{r){r 2 6ij - xiXj)dV (2.106) 

where p is the body’s mass density. An important, yet obvious prop¬ 
erty of the inertia tensor, is that it is symmetric. From this, it follows 
that all its eigenvalues are real with orthogonal eigenvectors (if the 
the eigenvalues are distinct the eigenvectors are orthogonal, otherwise 
we can choose an orthogonal basis for the degenerate eigenspace) and 
that by choosing a new coordinate system called the principal axes 
given by the three orthogonal eigenvectors, we can diagonalize the 
inertia matrix J, so that only the diagonal principal moments of 
inertia remain, with all the off-diagonal products of inertia van¬ 
ishing. For this reason, we will almost always choose the local body 
coordinate frame to align with the principal axes of inertia. For a 
body with finite volume, one can shown that in fact J is positive- 
definite, as well as symmetric. This ensures that all eigenvalues are 
positive and that our final equations of motion will have a unique 
solution. 


Examples of the Inertia Tensor: 

1. As a first example, consider a spherical mass of uniform density. Then 
from symmetry, the only non-vanishing elements of J are 

Jl.l = Jf 3 
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J 2,2 ” «//3 

(2.107) 


CO 

II 

CO 


where 

r 



J = 2 J p{r)r 2 dV (2.108) 


2 . A second example is fche symmetrical top , where the body has a ro¬ 
tational symmetry about some axis. Since we are free to choose the 
coordinates so that this axis is the z-axis, we have 

Jl t i = 

= (1/2) J p(f)(r J + z 2 )dV (2.109) 

Js,3 = P{r){^ 2 + y 2 )dV 

A large collection of examples can be found in (45). 

2.4.3.1 Newton-Euler Equations 

The motion of a rigid body can be decomposed into two parts: The 
translational motion of its center-of-mass, treated as a particle with 
the total body mass, and purely rotational motion about the center- 
of-mass. (A proof of this can be found in any mechanics text, e.g. 

130].) 

If r is the body’s center-of-mass position in some inertia! frame 
and m is its total mass, the equation of motion for r is just Newton’s 

m~r = f (2.110) 

where / are total external forces acting on the body. Suppose L is 
the body’s angular momentum. Then we know that 

L = Jw (2.111) 


where I is the moment of inertia of the body and w is the body’s 
angular velocity. This relation is true in any coordinate frame, in 
particular, in a frame fixed with the body. But Newton’s equation 
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for the rate of change for angular momentum expressed with vector 
components in the inertial frame gives 

L = n (2.112) 

(n the applied torque) when combined with our rate of change formula 

L — Lbody "1" wL (2.113) 


expressed in body-fixed coordinates gives Euler’s Equation 

J'u + wW = n' (2-114) 

since the inertia matrix J 1 is constant in the body frame. The diffi¬ 
culty with using this equation is that although J' may simple in the 
body-fixed frame (if we choose principal axes), the angular velocity 
vector is more complex. Often it is most convenient to express the 
external torques in terms of the local body frame. Equation 2.114 in 
component form is 

JjWj uj 2 lo 3 [J 2 *^ 3 ) — n \ 

= n 2 (2.115) 

j< .< ' > , j> x'v > 

^3^3 “ ^ 2 ) ~~ n 3 

where the vectors are expressed in the principal axis frame of the 

body. 

We can transform Equation 2.114 to global coordinates by making 
use of the cosine matrix A between the local body frame and the global 
frame. In this case, A is itself a function of the position of the body 
parameters, and A is not zero, 


JCJ - JAA t w + ujJCo = n (2.116) 


Example of torque-free motion of a rigid body: If there are no torques 
applied to the body, the angular momentum is conserved and the equations 
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for the angular velocity in the body-fixed principal axis frame are Equations 
2.115 

J 2<*>2 — <*>3Wl(>/3 — 7i) (2.117) 

JzOJS = U>i<a>2{Ji — J 2 ) 

We know that rotational kinetic energy and total angular momentum are 
conserved and so they are “integrals of the motion”. Using this fact, these 
equations can be explicitly integrated in terms of elliptic functions. There 
is also an interesting geometric interpretation of the motion due to Poinsot 
that can be found in [30], where the body’s motion can be reduced to the 
rolling of its inertia ellipsoid on the invariable plane whose normal vector is 
the fixed angular momentum L. During this rolling, the path traced out on 
the ellipsoid is called the polhode, while the corresponding path on the plane 
is called the herpolhode. Equations 2.117 are equivalent to the wonderful 
statement that “the polhode rolls without slipping on the herpolhode lying 
in the invariable plane”. 

If we assume a symmetric top, where 

J 1 = J 2 = J\2 # J 3 , 

then Equations 2.115 become 

Jl2Wl = <^2^3(7^12 — J 3 ) 

Jl2<*>2 — 0J3Wl(«/3 — 7 12) 

J3UI3 = 0 

Then 0 J 3 = constant and the first two equations become 

Jz ~ Jn 

=---W 3 W 2 

J\2 

Jz - 7 12 

o>2 = -:-W3 W 1 

7l2 

If we define f2 = then we have 

W] + HU2 — 0 
d>2 — fM = 0 


(2.118) 


(2.119) 


which has the well known solution 


( 2 . 120 ) 







Figure 2.12: Symmetric top. 


The vector w then has constant length and the angular velocity vector u> in 
the body fixed frame rotates in the , y'-plane with angular velocity fi. In 
the global frame, this corresponds to precession of the z'-axis. See Figure 
2.12 Equations 2.110 and 2.114 can be combined into a single matrix 

equation 

MY = F 

(2.121) 

where 

m 0 0 0 \ 

.. 0 m 0 0 

M ~ 0 0 m 0 

^ 0 0 0 J' j 

(2.122) 

is symmetric and positive-definite, 


Mi) 

(2.123) 

F=( 1 ) 

if- Ci'J’S' ,/ 

(2.124) 
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Figure 2.13: Two degree-of-freedom car model. 

The combined equations 2.121 for the translational and rotational 
motion of a rigid body is often called the Newton-Euler Equation. 
It is a first order equation in the translational and angular velocities 
v and c o. To calculate the body position and orientation, we add the 
coupled first-order system for translation 

r — v 

mv = f (2.125) 


Example of 2 degree-of-freedom car model: Consider the car coordinates 
shown in Figure 2.13. We assume it is traveling forward at a constant velocity 
U and with small angular deviations. Let x,y be the chassis c.g. coordinates. 
The front axle is a distance a in front of the c.g. and the rear axle a distance 
b to the rear. The tires are attached at a distance t from the centerline. The 
coordinates of the front wheels are 

xi t 2 = x + acos(t/>) ± t sin(i^) 

= x + a ± tip 

yi,2 = y + a sin(V>) =F 1 cos W 
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= x + aip T t 

The sideslip angle a, which is the angle between the tire hub velocity 
vector and the tire heading, is a = y/x-rp. Lateral complicated functions of 
normal force and slip angle. We assume that for small slip angles, the lateral 
forces are linear functions of the slip angle, with c/, c r the initial slope of the 
lateral force/sideslip curve at a fixed normal force. Since we have assumed 
that x = U, the slip angles for the front tires are atij = at/ = (y+aip)/U-ip. 
A similar analysis for the rear slip angles gives 03,4 = o r = (y - bip)/U — ip. 
Thus the Newton-Euler equations for the car motion are 

my = CfCtf + c r a r + Y ext 
hip = acfOtf - bc r a r + N ext 

where Y ext is the external force in the y-direction and JV ext is the external 55- 
moment. The combination of these equations with the equations for the slip 
angles provides a complete set for the determination of the vehicles motion. 


We have to choose a parameterization of orientation to unwind 
Euler’s equation into first order form. If we choose Euler parameters, 
then we must express u> in terms of the e’o, e, e 0 , t. From Section 2.3.2.4, 
we know that the orientation matrix M can be written in terms of 3x4 
matrices E and G involving the Euler parameters. Thus we could re¬ 
write Euler’s equation explicitly in terms of Euler parameters using 

Q = 2E ^ j (2.126) 

to get a second-order equation in the Euler parameters. This would 
require the addition of the normalization constraint between the Euler 
parameters. 

Instead of using the Newton-Euler equations directly, we will use 
the principle of virtual work which is much more adaptable to linked 
multi-body dynamics. We will separate the integration of the angular 
accelerations from the integration of the time derivatives of the Euler 
parameters. We use 


(f)=(1/2 jgv 


( 2 . 127 ) 
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to determine the Euler parameter time derivatives. We combine this 
Equation 2.127 with the Euler equation to find the body orientation. 

2.4.3.2 Virtual Work Equation for a Rigid Body 

A virtual work principle is true for a rigid body if we add to the 
translational work of moving the center-of-mass the work of rotating 
the body about the center-of-mass by a virtual rotation angle Sir 
against any external torques n ext . Because we consider only virtual 
rotations, they are infinitesimal and so are themselves vectors. All 
virtual displacements used in the virtual work principle must be in the 
same sense we have used before kinematically admissible. Also, just as 
in the particle case, we can ignore any internal reaction torques since 
we will assume they do no work. From d’Alemberts principle, we must 

also add the inertial terms from the rotational accelerations if r ( J w + 
wJV. The total additional rotational virtual work term we must add 

. I 

is Sn T (j'Cj +Cj j'— n' ext ). Notice that we use global coordinates for 
the center-of-mass motion given by Newton’s terms and body fixed 
coordinates for the Euler terms. To save writing the we will abuse 
the notion (and the reader) and drop the explicit reference to the 
body-fixed frame for any rotational term in the virtual work equation. 
This includes the virtual rotation H . It will be understood that these 
terms are always given in the fixed body frame. 

This gives the virtual work equation for a rigid body 

SZ T {Mf - F ext ) = 0 (2.128) 

where the virtual displacement consists of a kinematically admissible 
virtual translation Sr and virtual rotation Sir 

6Z = ( % ) (2.129) 

where again, F ext are the external forces and torques acting on the 
body. As in section 2.2.3, we can ignore workless internal constraint 
forces and torques. 


Example of virtual work with constraint condition: Suppose we have a 
wheel rolling in one dimension without slipping. Let 0 be the rotation of the 
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wheel from the initial vertical axis. Let x be the hub position of the wheel. 
The condition for roll-without-slip is dx + rdO = 0, where r is the tire radius. 
If the tire has mass m and moment of inertia about its polar axis I, then the 
virtual work equation with Lagrange multiplier due to the constraint is 

(mi - f ext + A)$x + (10 - n ext + r\)60 = 0 (2.130) 

where the virtual displacements for position and angle are independent since 
we have introduced the multiplier. Since the variations are independent the 
equations decouple to become 

mi = /"‘-A (2.131) 

10 = n ext + rA 


If we differentiate the constraint equation twice and substitute into the equa¬ 
tions of motion to eliminate 0, we can solve for A to get 


mrn ext + If ext 
mr 2 + I 


(2.132) 


If the tire is accelerating due either an applied external force or torque there 
is a non-zero multiplier A which is the frictional reaction force needed to 
keep the tire rolling without slipping. 


A collection of rigid bodies will have for their virtual work equation 


Y,&Zf{MiX i -F: xt ) = Q 

i 

where t refers to the body number and 



(2.133) 


(2.134) 


As with collections of particles, if constraints exist between the 
rigid bodies, the virtual displacements Z{ are not independent, and 
so the individual body Newton-Euler equations do not decouple. We 
must, as with the particle case, introduce Lagrange multipliers into 
the equation in order to assume the 8Z\ are independent. 
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2.4.3.3 Recipe for Multi-body Equations of Motion 
Following ([37]), we give a recipe for multi-body equations of motion: 

1 . Define the generalized coordinates q and determine whether they 
are dependent or independent. 

2 . Write each body’s Cartesian c.g. position r, = ^( 9 ) and ori¬ 
entation matrix Ai = Ai(q) as a function of the generalized 
coordinates. 

3 . Write the global body Cartesian velocities r, = and body 

| • * 
local angular velocity = Bi{q)q using vector angular velocity 

relations orw = A T A. 

4 . Write virtual displacements and rotations as 6 r,- = U^Sq and 
6G}\ = Bi{q)6q 

5 . Write the global body Cartesian accelerations r,- = r^q and 

. » - . ^ 

body local angular velocity = B,q + B,q using vector angular 
velocity relations or Co = A T A. 

6 . Derive applied forces f- xt and torques n txi . 

7. Calculate masses m, and local frame inertia matrices J,. 


8 . Substitute steps 3-5 into the virtual work equation 


nb 


JLitrf[rmfi - f- xt ] 4- 6nT\J.Q i + Co J,w - n/ 1 ']} = 0 
1=1 

and rewrite as (S^fM {q)q — Q‘ xt {q, q) — R ext {q,<l) ] = 0- Here we 
have inserted the “”’ 8 , denoting the body-fixed frame, back into 
the equation for clarity. 

9. If the q are independent, 6q is arbitrary and M (q)q = Q ext {q,q) + 
R' xt (q,~q). 

10 . If the q are dependent, derive the constraint equation $( 9 ) = 
0 , the Jacobian $$•, and the constraint acceleration equation 
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11. The constrained equations of motion are 

(« 

The solution of these equations is non-trivial because of the ex¬ 
istence of the holonomic constraint $ = 0 which must also be satis¬ 
fied. These equations are not differential equations, but differential- 
algebraic equations. In addition, initial conditions must be chosen to 
be consistent with the constraint equations $ = 0 and $ = 0. 

We will see in the next sections how to implement these equations 
so that the algorithm can be efficiently used on a parallel computer. 
The trick is to decompose the vehicle kinematics in such a way that a 
recursive elimination of the outboard relative joint generalized coordi¬ 
nates can be accomplished. This requires the elimination of all closed 
loops in the kinematic connections and the substitution of addition 
constraints to the equations of motion. 
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Chapter 3 

Real-Time Dynamics 
Algorithms 


3.1 Previous Dynamics Modeling Tech¬ 
niques 

Computational dynamics is a current very active area of research [34]. 
Previous work in obtaining equations of motion for a complex system 
can be generally divided into two categories depending on choice of 
system coordinates [35]: In the first, no attempt to choose a minimal 
set of coordinates is made, opting instead for a maximal set of Carte¬ 
sian coordinates for each body [54]; in the second class, a minimal, or 
at least small, set of relative coordinates is used recursively to identify 
a body’s position [38]. The first method has the advantages of gener¬ 
ality and relatively simple equations - and the disadvantage of rarely 
being optimal for a given mechanism or vehicle. Commercial codes, 
such as DADS [68] or ADAMS [18], generate Cartesian coordinates 
for each body to which a maximal set of constraint conditions must 
be applied. This is similar to finite element methods in structural 
mechanics where a very large, but sparse, equation set results. In 
the second method, a small set of relative joint coordinates is used, 
and fewer constraints must be applied, but the equations are often 
extremely complex. A minimal coordinate description is very effec¬ 
tive in describing spacecraft and robots, since these systems naturally 


75 



76 


CHAPTER 3. REAL-TIME DYNAMICS ALGORITHMS 


form tree graphs. But a problem occurs when closed loops appear in 
the system graph, since the relative joint coordinates are then sub¬ 
ject to additional constraints from the closing of the loop, and so are 
not independent. Wittenburg [70] considered the idea of cutting the 
closed loop to form a tree graph, and adding Lagrange multipliers to 
the equations of motion to satisfy the now missing constraints. Once 
the system graph is a tree, recursive methods developed by Hooker 
and Margulies [38] for spacecraft, and Armstrong [4] and Luh, et. 
al. [46], for robotics, can be used to recover body accelerations. Bae 
[5], using these ideas in his thesis, derived an algorithm combining 
cut joints with recursion. In addition, he discussed issues relating to 
the natural parallelization of the resulting algorithm. Essentially, we 
have implemented Bae’s work, adding additional techniques for sta¬ 
ble integration and pipeline delay reduction to the parallel recursive 
algorithm. 


3.2 Derivation of Recursive Dynamics Equa 
tions for Vehicle Model 

3.2.1 Rigid Body Kinematics and Euler Parame¬ 
ters 

For convenience, we summarize material from Chapter 2. As we have 
seen, a free rigid body in space has six independent degrees of free¬ 
dom: Three translational for its center of mass (CM) and three ro¬ 
tational for orientation about its CM. Traditionally, Euler angles are 
used to describe a body’s orientation. But such a parameter choice 
has singular configurations where its inverse transformation fails to 
exist. To avoid this potential difficulty, we use four Euler parameters 
c 0 , Ci, e 2 , e 3 , describing the angular orientation of the body fixed coor¬ 
dinate frame with respect to the global coordinate frame, satisfying 
a single normalization condition e* -I- e\ 4- e\ + e| = 1. The Euler pa¬ 
rameters (which are unit Quaternions) are related to space rotations 
in the following way: If the rotation has unit axis vector f} and angle 
9 about the axis, then e 0 = cos(0/2) and e = (ci,e 2 ,e 3 ) T = »7sin(0/2). 
Euler parameters describe the unit 3-sphere in Euclidean 4-space, and 
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Figure 3.1: Collection of rigid linked together to form a close loop. 

so clearly have no singular configuration. For this reason, they are 
a convenient choice of parameterization when arbitrary vehicle ori¬ 
entations are possible. We combine the Euler parameters as a four 
component vector p = (eo,ei,e 2 ,e 3 ) 7 ’, but we must remember that 
rotations themselves do not commute and so are not vectors. Simple 
algebraic formulas allow the expression of direction cosine matrices, 
virtual rotations, and angular velocities in terms of Euler parameters. 

3.2.2 Kinematic Relations - The Vehicle Graph 

Figure 3.1 shows a multi-body system of rigid bodies linked together 
by revolute joints which allow only relative rotation about the joint 
axis between bodies. The bodies are not independent in their motion, 
but are constrained to move so that the holonomic relations between 
the coordinates of each linked body are satisfied. 

We can abstract this picture and describe the same multi-body 
system in terms of a graph, where the vertices (nodes) represent rigid 
bodies and the edges joints. This is shown in Figure 3.2. A directed 
graph is a graph where an ordering between vertices is given, so each 
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\ 1 



Figure 3.2: Graph associated to multi-body system. 


edge has a direction from the ordering of its endpoints. It is possible to 
be very abstract and prove many theorems about graphs ((61]) useful 
for automatic generation of code based upon the graph structure of 
the multi-body system. Here we will be naive and not worry about 
graph-theoretic formalities. Additional constraint conditions between 
bodies not from joints are indicated by a dotted line. Figure 3.6 in the 
next chapter shows the graph for the linked rigid bodies making up 
a modern sports sedan. The chassis is the base body or first body in 
the directed graph. The constraints imposed by the steering linkage 
are shown as dotted lines. The linkages themselves are not treated 
as distinct bodies, but as distance constraints. It is clear that any 
set of linked rigid bodies can be described in this way, since we make 
no assumption that the graph have any embedding properties in the 
plane. 

We can now order the bodies by numbering from some starting 
base body, and derive the kinematic relations between any two linked 
bodies. Suppose we consider the i th and j th bodies linked together 
as shown in Figure 3.3. Here we have indicated several coordinate 
systems: The global frame; body i centroidal frame; body j centroidal 
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frame; the joint frame ij between body i and j; the joint frame ji 
between body j and i. J is the joint vector between the ij and ji 
joint frames. The various orientation matrices Ai , Aj , , C fJ -, Cjj 

relating the frames are indicated. The global centroidal coordinates 
of body i are r* and the joint attachment vectors in the local body 
frame are s { j. The Ai are use to transform local vectors like the 
to global coordinates = A t s|j. 

We assume that the joint between the two bodies can be fully 
described by a set of joint coordinates qij. For example, in a revo¬ 
lute joint, the angle about the joint axis vector; for a spherical joint, 
the Euler angles or parameters between the ij and ji joint coordi¬ 
nate frames; for a translational joint, the length extension along the 
translational axis definition vector. 

The following relations are easy to deduce: 


Aj 

= AiCijA" jCji 

(3.1) 


-4 

= ri + Si t j + dij — Sj t i 

(3.2) 

Uj 

= Wi + UJiJ 

(3.3) 
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Tj — 7"t -S, y + djj Sji , (^-4 ) 

where the of,- are the global angular velocities and oj,j is the global 
joint angular velocity between the ij and ji joint frames. A linear 
relation between and the relative joint coordinate velocities q i f - 
exists. We write it as 

u i,j = 

where Hij must be calculated for each joint type. Then 

u}j = ujf -+■ Hijqij. 

Since s ti y is a vector fixed on body i, 

5|,y = UfiSij. 


For body j, 




= - 5 j.« w r 


In terms of the relative joint velocities, 


S},i — 

The derivative of the joint displacement vector dij is 

dij — &idij 4 * di jq t } (]i j ( 3 * 5 ) 


where is the Jacobian of the joint vector dij with respect to 

the relative joint coordinates qij. Putting all this together, we can 
write the relations between the velocities in terms of matrices 
and Bin as 



A similar formula hold for virtual changes in the coordinates 6fi 
and 6 7T,- 




+ BijiSqij 


(3.7) 
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The acceleration relations are 


(t)=Mt) +a "^ +D 

tj (3.8) 

where Dij contains the velocity terms 


Aj = A.,,i (|‘ ) + B u ,, lj 

(3.9) 

The matrices Bij t i and 


B .. = ( 0 ) 

1,1 ^ 0 0 / 

(3.10) 

are independent of joint type, while Bij j and Bijj depend on what 
the joint and its relative coordinates. Reference (5] has an extensive 
table of joint types and their Bi and f ?2 matrices. 


Example of revolute joint: Figure 3.4 shows to bodies linked by a revolute 
joint. Here dij = 0. Choose the joint axis vector h parallel to the z < ; -axis 
and let 0 be the rotation angle about this axis. Then 

iijj = &i + HijOij 


where 



(3.11) 

(3.12) 


The joint transformation matrix A"j is just the usual rotation matrix about 
the local ij joint frame by the angle 0. Also, 


and 


d . _ ( hihj 

1 hi 


•y,i*u + h<hi 

hi 


(3.13) 


B i,i, 2 = 


(3.14) 



82 


CHAPTER 3. REAL-TIME DYNAMICS ALGORITHMS 



Figure 3.4: Two bodies linked by a revolute joint. 


Since 


h-ij — 


D _ ( hjj&jSj'i + Sj'iUihi'j \ 
Si " - { ) 


(3.15) 


3.2.3 Virtual Work Equations 

The virtual work principle states that for a system of N bodies, the 
sum of the virtual work (work of a virtual displacement 6q against 
a force Q) due to inertial and external forces under virtual kinemat¬ 
ically admissible displacements (virtual displacements which satisfy 
the constraints), added to the the work of the internal forces, must 
vanish in dynamic equilibrium. For a rigid body we ignore this last 
term, because the internal forces can do no work. We treat springs 
and dampers connecting rigid bodies as external force elements. 

For a rigid body, we must consider the virtual displacements from 
both translation and rotation. We follow the notation of [37). If the 






83 


position and orientation are given in the global inertial frame by 
and pi, we write the i-th body’s global velocity as the vector Yi 


Yi 



(3.16) 


where its rotational velocity is w,-; the body’s acceleration and virtual 
displacement-rotation are 


Mi:)’ 

(3.17) 

«>) 

II 

•* 

(3.18) 


where the virtual rotation Sir is a true 3-vector. If the system of bodies 
must satisfy a set of constraint equations $(ri,p,) = 0, then for every 
set of kinematically admissible virtual displacements SZ, satisfying 
orthogonality with the rows of the constraint Jacobian where 

this Jacobian notation represents a matrix of partial derivatives of 

each component of $ with respect to each coordinate, expressing the 
linear variation of the constraints in terms of the virtual displacements 

AT 

EV^'= 0’ ( 319 ) 

i=l 

the virtual work must vanish: 

f^6ZT{M i 9 i -Qi)= 0, (3.20) 

i=i 


with the symmetric and positive-definite generalized mass matrix M< 
and force vector Qi given as 



mi 

0 

0 

0 

\ 


0 

m, 

0 

0 



0 

0 

m, 

0 



0 

0 

0 

Ji 

/ 


(3.21) 
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and 

Qi = ( h _ . . ) (3.22) 

\ n,- -Wj x JjW J 

contain the body mass m,-, inertia Jj, external forces and torques 
Hi. Geometrically, kinematic admissibility requires every virtual dis¬ 
placement be in the tangent hyperplane to the configuration manifold. 

In our case, the virtual displacements 6Zi are not independent 
when subject to the constraint conditions Equation 3.19, and so we 
cannot assume each of their coefficients must vanish. If we add the 
virtual work of these constraint forces, by means of Lagrange multi¬ 
pliers, then the virtual displacements can be assumed independent. 
So we replace Equation 3.20 by 

- Q.) + £ «#*£*%•.»)=<>, (3.23) 

i=l (j t k)—cut joint 

where the are the Lagrange multipliers from the cut joint be¬ 

tween body j and k. 

3.2.4 Recovery of Accelerations from the Recur¬ 
sive Equations 

We demonstrate (following Bae’s thesis (5], [6], (7], [8]) how the re¬ 
cursive equations for the accelerations are derived by considering the 
single closed loop graph shown in Figure 3.5. Body 1 is the base body, 
and we cut the loop between bodies 1 and 1+1. Equation 3.23 is then 

£ -<?<) + £ Zf (Mi?i -4 + «2;' +,) %,, + „) = o. 

»=i.«yu+i t=i 

(3-24) 

If we now substitute our kinematic relations relating acceleration and 
virtual displacement of body 1 in terms of those of body 1-1 and the 
relative joint coordinates between them, we get 


-Qi)+ 


1 = 1 
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Figure 3.5: Closed graph for derivation of recursive elimination equa¬ 
tions. 


SZf (Bj+ M,/?(,_!), - Qi) + ^g; I+1)T A (I . l+ i ) ))-h (3.25) 

- Qi + * ( J;‘ + 1 )T A (m+ i)))+ 

£ - Q,-) + - Q,+i + *S I ^° r A ( i 1 , + i,) = 0. 

*= 1+2 

Since the virtual relative joint displacements Sq^-iy are arbitrary (af¬ 
ter making the cut), its coefficient must vanish and we can solve for 
the relative joint acceleration 

(3.26) 

If we substitute the relative acceleration back into Equation 3.26, we 
have eliminated all body 1 accelerations. Continuing in this manner 
backwards along each branch gives an equation for the base body 

acceleration Y i and Lagrange multipliers A 

6 % [{Mi + K\ + K\)?i - (Qi + L\ + L\) + (*M, ir + = 0, 

(3.27) 
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where the superscripts 1 and 2 indicate the result of recursive elimi¬ 
nation along each branch to the base body for the accumulated mass 
(K), force (L), and Jacobian ($M) matrices. 

If the system has nc constraints, Equation 3.27 does not yet deter¬ 
mine the base body acceleration since we have only 6 equations and 
6 + nc unknown accelerations and Lagrange multipliers. To arrive at 
the extra equations needed, we differentiate the constraint equations 
twice to get 


$(,_!), = $ f( y, + $2 i+i yi +1 -7 

= 0 


(3.28) 


where 7 includes all other non-acceleration terms. We can now pro¬ 
ceed as before to eliminate recursively the intermediate bodies down 
to the base body, arriving at 

($M X + + ($L} + $L \)A - RHSl - RHSf -7 = 0 , (3.29) 

where again we have accumulated Jacobian mass ($M), force ($L), 
and mixed term (RHS) matrices along each branch. 

We can combine Equations 3.27 and 3.29 as 


(3.30) 


where we have collected terms from both branches. Equation 3.30 is 

then used to find the base body acceleration Y Once the base body 
acceleration is known, we recover each outboard body relative joint 
acceleration using Equation 3.26 and then the Cartesian acceleration 
using the kinematic formulas. 


( *1 

$Af x r \ 


\-f Lt ) 

\ *M x r 

) 

1 * j 

1 " ^ RHSx ) 
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3.3 Vehicle Dynamics Model 

3.3.1 Real-Time Dynamics for Vehicle Simulation 

Evans & Sutherland has devised efficient real-time dynamics algo¬ 
rithms applicable to arbitrary vehicle types. A description of how 
these algorithm work is fairly complex, but can be broken down into 
a number of areas. 

Non-real-time computer simulation of complicated mechanical link¬ 
ages has been performed for machines, robots, spacecraft, and ground 
vehicles. True real-time simulation is more difficult and hitherto done 
only for relatively simple mechanical systems, such as aircraft, when 
the governing equations can be approximated, as with linearized equa¬ 
tions of motion, or with specialized computer hardware. Engineer¬ 
ing simulation of modern vehicles, such as high performance sports 
cars, requires a multi-body non-linear model in order to provide a 
high fidelity of motion; software flexibility is required to deal with 
the general nature of vehicle engineering problems. A solution to 
real-time engineering simulation demands both cost-effective compu¬ 
tational capacity and the adaptability of general purpose computers 
and software. Simulator system integration and introduction into 
the real-time software of motion-base washout and vehicle subsystem 
models are facilitated by the use of standard operating systems and 
programming languages. 

Two significant approaches to real-time vehicle dynamics mod¬ 
eling, albeit at very different levels, are those implemented by the 
Swedish Road and Traffic Research Institute (VTI) [52] and Daimler- 
Benz [27] [33]. The VTI model considers a simple planar 3 degree-of- 
freedom car, but carefully considers tire and other forces. Daimler- 
Benz achieves real-time with a much more sophisticated, but still 
simplified model. They incorporate an implicit integration scheme 
with a quasi-Newton method to achieve stability in a dual processor 
system. 

We are primarily interested in the accurate prediction of vehicle 
dynamics, rather than kinematics or inverse dynamics [39]. The ap¬ 
proach we use decomposes the system into a set of rigid bodies (a 
rigid body requires distances between points within the body to be 
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Cut 



Figure 3.6: Vehicle graph with cut-joints indicated. 


fixed) connected by joints (which allow relative motion between the 
bodies), subject to various forces, moments, and constraints. The re¬ 
sulting linked system can be described by a directed graph [70] (see 
Figure 3.6), where the bodies form vertices and joints the edges join¬ 
ing the vertices. People and robot arms, usually without closed loops 
in their connections, form special tree graphs contractible to a single 
base vertex. A tree graph has a natural direction outward from the 
base vertex. In a tree graph, outboard body (directed out along a 
branch in the graph) Cartesian positions and velocities can be ex¬ 
pressed in terms of inboard body and relative joint coordinates and 
velocities, recursively providing the coordinates needed for determin¬ 
ing motion. We use the term recursive to mean either moving inward 
or outward along the tree graph, using data from previous bodies 
along the branch to determine the position and motion of the next. 
General graphs, with embedded closed loops, are required for more 
complicated systems, such as high-performance sports cars. Our kine¬ 
matic description of a vehicle allows use of graphs with closed loops. 
To derive equations of motion for our system, we first cut the closed 
loops, so that a recursive relative joint coordinate kinematic descrip- 
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tion can be employed. Next, D’Alembert’s virtual work principle [31] 
is utilized to derive the equations of motion, augmented by reaction 
forces from the closed loop joint constraints via Lagrange multipliers. 
The relative joint coordinate formulation, implemented in a recursive 
parallel algorithm, leads to greater efficiencies for real-time simula¬ 
tion. External forces are computed from the tires, drivetrain, wind 
resistance, brakes, springs, dampers, stabilizer bars, and suspension 
stops. Because our vehicle bodies are subject to kinematic constraints, 
by cutting the closed loops in the graph, we produce a set of coupled 
differential-algebraic equations (DAE’s), in which differential equa¬ 
tions of motion and algebraic equations of constraint must be simul¬ 
taneously solved. At all times, the numerical solution must closely 
satisfy the imposed kinematic constraints, so some form of constraint 
stabilization is required. In addition, stability must be preserved even 
with high frequency force inputs. 

In the rest of this section, we describe our recursive vehicle kine¬ 
matics and equations of motion. Next we treat the numerical in¬ 
tegration techniques required for our DAE’s, along with the impor¬ 
tant concept of constraint stabilization. Lastly, we mention computer 
hardware issues and present the timing results for our algorithm on 
various serial and parallel machines. 


3.3.2 Vehicle Description 

A typical modern high-performance sports sedan is shown schemati¬ 
cally in Figure 3.7. It consists of a chassis, the two MacPherson front 
suspension subsystems [15], the steering rack, and the two rear suspen¬ 
sion arms. The MacPherson front suspension subsystem is detailed 
in Figure 3.8. This suspension allows both jounce (spindle vertical 
motion) and steer (the rotation of the spindle) as independent mo¬ 
tions. Steering is done by motion of the tie-rods (which act as distance 
constraints) connecting the rack and spindle. 

The initial step is to break the vehicle into its independent struc¬ 
tural components, each of which can be approximated as a rigid 
body. These will include everything except compliant members, such 
as shock absorbers, bushings, and springs, which have important in¬ 
ternal degrees of freedom. These compliant elements can themselves 



90 


CHAPTER 3. REAL-TIME DYNAMICS ALGORITHMS 



Figure 3.7: Schematic diagram of modern high-performance sports 
sedan. 



Figure 3.8: Macpherson front suspension. 
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be modeled using deformation mode analysis, with the coefficient of 
the mode introduced as another generalized degree-of-freedom. For 
example, a model might contain 10 independent bodies: Chassis, rack, 
two Swiss arms, two struts, two spindles, and two rear trailing arms. 
Shock absorbers and springs are attached to the front and rear sus¬ 
pensions. Torsion bars are included to provide the proper vehicle roll 
stability. 

Each rigid body needs six coordinates to describe its motion: 
Three for its center of mass and three to describe its angular orien¬ 
tation. Because the rigid bodies are connected to one-another, they 
are constrained in their motion. For example, the Macpherson front 
strut assembly has only two motions: The vertical “jounce” and ro¬ 
tational steering. When we account for all the constraints from the 
joints, we have just ten independent degrees of freedom left. The joint 
constraints effectively convert the dependent individual body coordi¬ 
nates into a smaller set of independent relative joint coordinates. For 
example, the rear trailing wheel only needs the relative angle between 
it and the chassis to describe its position. The model takes as its basic 
variables the chassis Cartesian coordinates and relative joint coordi¬ 
nates. These, along with the associated velocities, are the coordinates 
that must be propagated through time by the numerical integrator. 

3.3.3 Real-time Algorithm 

3.3.3.1 Decomposition in terms of a Vehicle Graph 

The total vehicle is describable by a graph structure: Each rigid body 
is a node, and the nodes are connected by directed edges which rep¬ 
resent the joints. In general, closed loops will exist within the graph. 
For example, in the front suspension, a closed loop occurs among the 
chassis, Swiss arm, spindle, and strut. These closed loops will cause 
problems for the relative coordinates, so they are removed by apply¬ 
ing a “cut-joint” constraint in place of an actual joint. This produces 
a tree-structure graph, without closed loops, but adds additional con¬ 
straint equations that must ultimately be satisfied by the coordinates. 

The corresponding vehicle graph is shown in Figure 3.6. There are 
four closed loops from the front suspension and steering linkages. If 
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we cut the links between rack and spindle and lower control arm and 
spindle on both sides, we produce the tree graph shown in Figure 3.9. 


3.3.3.2 Vehicle Coordinates 

In Figure 3.9, the chassis is the base body, whose coordinates are 
the three Cartesian degrees of freedom of the center-of-mass in the 
inertial frame r\ and a vector 7?! of four Euler parameters describing 
its orientation. The lower control arm position is determined by the 
chassis position and specification of the revolute joint angle (the rear 
suspension arms have a similar description). The strut’s position is 
given by the chassis and spherical joint orientation (again using Euler 
parameters). The spindle’s position can be found from the strut’s 
position and the translational joint distance. The position of the 
rack relative to its translational joint with the chassis is given by the 
driver’s steering input. Table 3.1 presents the system’s coordinates. 
Notice that we have not yet “closed” any of the four loops in Figure 
3.9 by adding constraints terms to the virtual work equations. 
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Table 3.1: Model coordinates. 

The four tires are not treated as individual bodies, but are given 
only a single rotational degree of freedom and inertia. In our model, 
there are 10 rigid bodies, 4 tires, and 26 body and relative joint coor¬ 
dinates (not all independent) used in the kinematic description of the 
vehicle tree graph. The constraints from the closed suspension and 
steering loops must be added to these 26 coordinates. Not treating 
the tire as a full body is a reasonable approximation, since tire inertial 
forces can be separately calculated and added as additional external 
coupling torques to the equations of motion. Compliance in any of 
the bodies is not considered (cf. remark at the end of the dynamics 
section). 


3.3.3.3 Holonomic Vehicle Constraints 

The coordinates used in the vehicle description are subject to holo¬ 
nomic constraints [30] $(J?<) = 0 (constraints in precisely this func¬ 
tional form, which can be used to eliminate dependent coordinates) 
from the closed loops in the graph. The MacPherson front suspen¬ 
sion requires the spindle to meet the lower control arm at a spher¬ 
ical joint as shown in Figure 3.8. Counting both the right and left 
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sides, this gives 6 independent constraints. The rack-spindle distance 
give 2 more independent constraints, while the normalization of the 
three sets of Euler parameters provides 3 more. The driver steering 
commands produce 1 more constraint on rack motion. Since each of 
these constraints is holonomic, and they are all independent, we have 
26 — (6 -I- 2 + 3 -|-1) = 14 independent degrees of freedom in our vehicle 
(exclusive of the drivetrain and other subsystems). 

Non-holonomic constraints (non-integrable relations involving the 
coordinates) are encountered in our model only in the kinematics of 
tire roll-without-slip on a surface. The resulting tire constraint forces 
are treated as external forces on the spindle and tire. 

3.3.3.4 Recovery of Outboard Positions, Velocities, Accel¬ 
erations 

We work with a set of 26 dependent coordinates describing the tree 
graph structure and 12 constraint equations from the cut joints. The 
velocity and acceleration formulas involve joint matrices B i and Bz 
giving the linear relation between outboard body velocity and inboard 
body and relative joint velocity. We write the Cartesian velocity of 

outboard body i+1 as Zj +1 = B x Zi + 5 2 7( 1+1 j, where for body i, its 
global Cartesian coordinates are Z, = (rf nf ] T . . The corresponding 
acceleration formula is Zj+i = B X Z{ 4- B29(< + i). + where the 

velocity coupling terms Z)(, +J )i contain the coordinate velocities and 
derivatives of the joint matrices B\ and Bz- The B{ and D matrices 
must be individually constructed for each specific joint type. 

These kinematic formulas are used to express the global Cartesian 
position, velocity, and acceleration of an out-board body in terms of 
an inboard body and relative joint position, velocity, and acceleration. 
All body coordinates, velocities, and accelerations are determined re¬ 
cursively moving outward in the graph from the initial base body. 

3.3.3.5 Recovery of Coordinate Accelerations 

We generate our dynamical equations using D’Alembert’s principle of 
virtual work. This is done in a way that generates maximum paral¬ 
lelism in the algorithm. By using relative joint coordinates, rather 



than Cartesian coordinates, the unknown joint constraint forces are 
automatically eliminated from our equations. Using virtual work has 
the advantage that alternative variational methods require the addi¬ 
tional step of forming either Lagrangians or Hamiltonians, and so are 
cumbersome for complex systems. Added difficulties with variational 
formulations occur when non-conservative forces (such as friction) are 
present. With virtual work, we can consider both conservative and 
non-conservative forces in a straightforward manner, while the cut 
joint constraints are easily incorporated using Lagrange multipliers. 

The result of the equations of motion is a parallel algorithm that 
allows the inward recursive elimination of the relative joint acceler¬ 
ations, and accumulation of terms into a single matrix equation for 
the base body accelerations and Lagrange multipliers. The remaining 
coordinate accelerations are then recursively recovered by moving out¬ 
ward in the graph. The organization of the algorithm closely follows 
the system tree graph, and so exhibits the inherent parallelism of each 
independent branch. Accelerations are numerically integrated to de¬ 
termine the new positions and velocities at the next frame time. The 
recursive algorithm is fully non-linear in that velocity coupling terms 
D are not assumed negligible. The assumption of slow rates of angular 
change, so coupling from the Euler terms w x J u7 is small, fails during 
high speed maneuvering of vehicles. The inclusion of these non-linear 
terms is critical for simulating the effects of anti-lock braking systems 
(ABS) [41] and four-wheel steering controls. 


3.3.3.6 External Forces 

We include various external forces such as tire forces, wind resistance, 
engine and brake torque, springs, dampers, stabilization bars, and 
internal reactions in the model. The most critical, and difficult, for 
accurate road vehicle simulation are the tire forces [20]. The major 
components of tire forces and moments are shown in Figure 3.10. 
Because of the complexity of interaction, the forces and moments are 
taken from measured data as a function of normal force, slip angle, and 
camber angle. Attention is paid to how the tire forces change during 
conversion from a no-slip to a slip condition in order to minimize 
longitudinal force discontinuities. 
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Z (Normal Fore* F z ) 



Figure 3.10: Tire forces and moments. 

Complex drivetrain models for both manual and automatic trans¬ 
missions have been developed. They are based on combining driv¬ 
etrain components such as engine, clutch, transmission, torque con¬ 
verter, and differential together with proper feedback (50). Each sep¬ 
arate component is modeled using table lookups of data provided for 
individual designs. Our current drivetrain models have a single de¬ 
gree of rotational freedom which is independently integrated to find 
rotational velocity and applied wheel torques. 

3.3.3.7 Integration of Accelerations 

We can symbolically write the result of the recursive algorithm deter¬ 
mining the accelerations as recovering the right-hand side f(t,Y) of 
the first-order equation 

? = f[t t Y), (3.31) 

where Y is the state vector of 26 chassis and joint coordinates and 
their velocities. This equation includes both the dynamics equations 
and the differential form of the algebraic cut joint constraints. The 
resulting DAE’s are stiff and awkward to solve (28). Equation 3.31 
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can be solved usinj standard techniques from numerical analysis to 
advance K(t) to Y(t 4- At), if combined with constraint stabiliza¬ 
tion. Because of real-time limitations, methods that require only a 
single updated evaluation of / for each time step can currently be 
used. This limits us to conditionally-stable explicit methods such as 
Adams-Bashforth and Adams-Moulton Predictor-Corrector schemes 
[43] without a final update. Another option is to use low-order Runge- 
Kutta methods. These have adequate accuracy, with enlarged stabil¬ 
ity regions. If explicit methods are used, great care must be taken in 
how the resulting equations are partitioned and solved to avoid insta¬ 
bilities resulting from high frequency force inputs. To improve total 
system bandwidth, after the system has been properly partitioned into 
slow (chassis) and fast (suspension) components, a hybrid integration 
scheme where the low frequency slow components are integrated with 
an explicit method and the high frequency fast components with an 
implicit method is advantageous. 

Guidance in partitioning the components can be gained by the 
use of an eigenvalue analysis of the Jacobian matrix for the differen¬ 
tial equation system. This linear analysis, though only approximate, 
will help reveal the potentially troublesome high-frequency oscilla¬ 
tory modes (those with large imaginary parts) and stiff, fast decaying 
modes (those with large negative real parts). When combined with 
a stability region analysis of the ODE method, an informed guess as 
to stable step size can be made. Since it is very difficult to calculate 
the system Jacobian analytically, a difference approximation about a 
reference state of interest is used. The zero eigenvalue, rigid body 
modes must be first removed from matrix before EISPACK or other 
standard eigenvalue software is used. 


3.3.3.8 Constraint Stabilization 

The numerical solution to Equation 3.31 must satisfy all cut joint 
constraints. These additional algebraic equations (forcing us to solve 
the differential equations on the curved configuration manifold, rather 
than in flat Euclidean space) require special techniques in correcting 
the solution produced by the integration method, while retaining ac¬ 
curacy and stability. For short times, since the differential form of the 
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constraints are incorporated into Equation 3.31, the solution will ap¬ 
proximately satisfy the cut joint constraints. For longer times, either 
Baumgarte’s method [9) of modifying the equations to make the con¬ 
figuration manifold (the manifold determined by the constraints on 
which the system trajectories move) a stable attractor (points close 
to the configuration manifold move nearer in time along their tra¬ 
jectories), or a coordinate partitioning method (68] of updating the 
dependent coordinates can be used. Baumgarte’s method essentially 
introduces a feedback loop to control error growth, while the partition¬ 
ing method uses Newton-Raphson iteration to make the dependent 
coordinates satisfy all constraint conditions. Either of these ensures 
that accumulated numerical error will not cause the solution to drift 
too far off the configuration manifold. 

3.3.4 Organization of Parallel Recursive Algorithm 

The efficiencies required by real-time demand a careful decomposition 
of the algorithm data flow and attention to mapping the algorithm to 
a particular computer architecture. You can’t fit a shoe that doesn’t 
fit; you can’t expect real-time performance with algorithms and ar¬ 
chitectures that don’t match. 

3.3.4.1 Algorithm Flow 

Figure 3.11 shows an implementation of the recursive algorithm with¬ 
out constraint stabilization. At the start of each time step, Cartesian 
positions and velocities of each body are recovered from the relative 
joint positions and velocities moving outward in the graph from the 
chassis base body. Each separate branch is computed concurrently 
with the others. The driver steering, brake and accelerator pedal 
position and other inputs are obtained from sensors and used to de¬ 
termine rack position and velocity, brake, and drivetrain parameters. 
External forces and moments are next calculated; tire, brake, spring, 
damper, stabilizer bar, aerodynamic, gravitational, and drivetrain are 
independently determined by special subsystem models. Next, the 
accelerations are recursively eliminated moving inward in the graph 
toward the base body. Again, these eliminations are performed con- 
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Figure 3.11: Organization of parallel tasking in recursive algorithm. 
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currently. When the base body Equation 3.30 has been assembled, it 
is solved for the chassis Cartesian accelerations using Gaussian elim¬ 
ination, and recovery of the outboard body Cartesian and relative 
joint accelerations is made moving outward in the graph. When this 
is completed, the right hand side of Equation 3.31 has been calcu¬ 
lated and the integration scheme advances the solution to the next 
time frame. 

3.3.4.2 Computer Architectures and Software 

The requirement of implementing a parallel algorithm with a small 
number (4-8) of concurrent tasks leads to considering a single bus 
or cross-bar architecture with a shared global memory. Care must 
be taken to avoid potential bus contention from simultaneous mem¬ 
ory access or message passing. Modern optimizing FORTRAN and C 
compilers running under POSIX-type operating systems provide soft¬ 
ware flexibility. Specialized simulation computers with non-standard 
languages suffer from insufficient adaptability, limited transportabil¬ 
ity, and narrow growth potential. Upgrading of a simulator computer 
system is simplified by writing software in standard languages. 

A solution to these problems is to choose a standard vector-concurrent 
mini-supercomputer, like the Alliant FX/8, running a UNIX operat¬ 
ing system. Figure 3.12 shows the architecture of the Alliant FX/8 
(| 2 |)- 

In the Alliant FX/8, up to 8 compute elements (CE) are connected 
by a 376MB/sec cross-bar to two large computational processor caches 
(CPC). Each CPC is a two-way interleaved 256KB cache for up to 4 
CEs. Cache coherency is maintained by a write-back cache protocol 1 
watching the main memory bus. The CPC are connected to a 4-way 
interleaved 64MB main memory by two independent 188 MB/sec 72 
bit (64 bits data and 8 bits SECDED) data buses, a 28 bit address 
bus, and a concurrency control bus. Each individual CE contains an 
integer, floating point, and vector register set. The floating point units 

‘A write-back cache is one where writes to memory are stored in cache and 
written to memory only when a rewritten item is removed from cache. In a write- 
through cache, writes to memory are written concurrently to both the cache and 
memory. 





CONCURRENCY 

CONTROL 
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are fully segmented, and can perform simultaneous scalar multiply- 
adds. Special vector operations are available that will chain standard 
vector calculations using adds and multiplies together. All floating 
point operations are in 64 bit, and the vector lengths are 32 64-bit 
vectors (CRAY has a length of 64 64-bit vectors). The peak “macho- 
flops”, with all 8 CEs, are 188 MFLOPS. A more realistic number 
would be 20-40 MFLOPS sustained. 

The operating system and all I/A operations are done with in¬ 
teractive processors (IP). These are essentially Motorola 68020 units 
with some possibility of doing floating point operations. These IPs 
access main memory through the IP caches (IPC), each of which share 
up to 3 IPs. Each IP also has a local memory to hold operating system 
kernel text pages. Direct memory access (DMA) is available to VME 
devices through special VME IPs. The design of the Alliant allows 
computation to occur concurrently with I/A operations. If the major¬ 
ity of the executing process is cache resident, I/A operations have no 
effect on system performance. This is important for a simulator, since 
a modest amount of (less than 100 bytes) output must be delivered 
without slowing the real-time computations. 

Concurrency control is done through the existence of concurrency 
registers on each CE, and global control of the concurrency is done 
by concurrency instructions communicating via the concurrency bus. 
Synchronization of parallel execution of loops does not require access 
to main memory by the use of spin-locks and semaphores, but is done 
directly in the concurrency registers. 

Very good optimizing vector-concurrent FORTRAN and C com¬ 
pilers are available for this class of architecture. They are limited in 
their detection of parallelism by the inability to fully trace all data 
dependencies, and so really only automatically produce COVI - con¬ 
current outer, vector inner - loop breaking. If the algorithm does not 
allow for its description in terms of single or multiple loops, the pro¬ 
grammer must provide explicit directions to enable the compiler to 
parallelize the code. If Execution of multiple concurrent threads is re¬ 
quired, a time consuming analysis of the algorithm data flow must be 
performed, resulting in a decomposition into serial groups of parallel 
tasks. The compiler can then proceed to automatically insert instruc¬ 
tions for the required data coherence and task synchronization. The 
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RELATIVE TIME TO FINISH ONE ALGORITHM STEP 
(ALLIANT FX/8 WITH 6 CE’S - 1 UNIT) 


Alliant FX/8 

(8 CE’s - parallel code) 

1 

Alliant FX/8 

(4 CE's - parallel code) 

1.25 

Alliant FX/8 

(2 CE's - parallel code) 

2.0 

Alliant FX/8 

(1 CE - serial code) 

3.2 

VAX 8800 

(serial code) 

4.5 

VAX 780 

(serial code) 

21.5 


Table 3.2: Timing results for recursive algorithm. 


goal of the algorithm analysis and optimization is to provide as many 
task groups of up to 8 course-grained (to reduce the need for mem¬ 
ory conflicts over shared data) concurrent tasks as possible, with the 
smallest possible number of single serial tasks separating the concur¬ 
rent task groups. Such a design will maximize the parallel speed-up, 
hopefully approaching one linear in the number of processors. Practi¬ 
cal algorithms rarely achieve linear speed-up. The recursive algorithm 
on the Alliant FX/8 achieves a speed-up of 3.2 with 8 CEs (cf. Table 
3.2). 


3.3.4.3 Algorithm Timing 

The recursive algorithm for the vehicle model described in section 

3.3.3 has been run on the Alliant FX/8 and VAX-780. 

Table 3.2 shows the results of timings taken on these computers, 
both serial and parallel, using optimized FORTRAN codes for each 
machine. 
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3.3.4.4 Remarks 
Two remarks: 

1. Notice that the way the algorithm proceeds, a correction to tire 
position occurs. We do not just assume a tire position, but 
a step is taken to correct it, and improve the reaction force 
calculations. 

2. The algorithm, being organized around a tree structure without 
closed loops, clearly lends itself to a parallel implementation. 
The left-right symmetry produces the greatest improvement in 
speed, with additional marginal improvements among each set 
of branches on each side. 



Chapter 4 

Physiology for Simulation 


4.1 Effective Simulation 

Effective simulation requires that the simulated environment approach 
that of the real world as closely as possible. The various senors of the 
human subject need to be stimulated in such a way that proper cor¬ 
relation between all the sensors is achieved. Consideration of human 
factors in the design of the visual and motion systems are critical for 
a believable environment. 

Achieving the desired environment requires an understanding of 
the human perception system and its interaction with the simulator. 
One needs to understand the physiology of the various sensory organs 
to be able to coordinate the stimuli given to the subject so that valid 
results can be obtained. A simulator is of little use if it elicits an 
improper response to the presented scenario. 

First, a review of linear systems theory will be presented as it 
pertains to modelling the humans sensory systems and motion control. 
Next, a discussion on the human perception systems will be presented 
followed by issues regarding the design of an effective simulator motion 
system. 


4.1.1 Linear Systems Theory 

In subsequent sections, linear systems theory will be used to describe 
several different systems. We want to provide a brief review of the 
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theory and then illustrate the concept with an example. 

The first step in describing a system is to derive a set of linear, 
constant coefficient, ordinary differential equations (ODE). This set 
of equations describe the dynamics, or behavior, of the system. 

The first tool that will be used is the Laplace transform. In this 
transform we replace derivatives with the Laplace variable, s. The 
main advantage is that this transform puts the ODE’s into the fre¬ 
quency domain in a manner similar to the Fourier transform and 
transforms them into algebraic equations. These can then be ma¬ 
nipulated with the ease of any algebraic equation. All the standard 
linear system tools, such as frequency response graphs, are available 
for system analysis and design. 

The next step is to define the system transfer function. The trans¬ 
fer function defines the input-output relationship of a linear, constant 
coefficient system. It is important to note that it is only defined in 
the Laplace domain. Once the system transfer function is derived, the 
large toolbox of linear systems analysis tools are available to describe 
and characterize the original system. 

If the reader wants a more in depth review of linear systems we 
recommend Dorf, [25] or Kuo, [42] for control systems. For digital 
filters we suggest Oppenheim, [53] and for general linear systems, 
Strang [66] is a good reference. 


Example of the Laplace transform method for the differential equations 
of a spring, mass, damper system: 

The equation of motion for a generic spring, mass, damper system is: 

+c^ + Ky = r (t) (4.1) 

Where M is the mass, c is the damping, K is the spring constant, y is the 
displacement of the mass, and r(f) is the forcing function. 

Rearranging into the more familiar form with the natural frequency, 
u>„ = \ the damping ratio, f = ~, and the critical damping, c c = 2M«„: 


d 2 y 

dt 2 


+ ~ + w nV = 


(4.2) 
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The Laplace transform of 4.2 is: 

s J T(s) + 2fw n sV(s) + w*V(s) = (4.3) 


Where T(s) is now the transformed displacement, $ is the Laplace variable, 
and we have assumed that the initial conditions are zero. The transfer 
function is found by forming the ratio of the terms multiplying Y (s) and the 
terms multiplying i?(s). 


Y{») _ G(m\ _ 

R(s) ' K(a 2 + 2fw n s + w*) 


(4-4) 


Not all systems can be adequately described by linear systems the¬ 
ory. Even systems with mild non-linear behavior can display behavior 
that is not a small perturbation of some linearized model. 

4.1.2 Human Perception Systems 

The human has many sensors to allow us to define and interact with 
our environment. Proper stimulus of these sensors is the key to creat¬ 
ing a realistic environment. A proper understanding of the function 
and modeling of these systems is essential to their proper stimulus. 
The three basic groups of sensors that will be discussed are: 

1. The eyes. 

2. The vestibular system or middle ear. 

3. The proprioceptive sensors. 

4.1.2.1 The Eyes 

The eyes are perhaps the most complicated system we have for per¬ 
ception. In the simulator, the eyes are used as a primary source of 
information about the environment. A complete exposition of the eye 
and its’s functions is well beyond the scope of this paper. The reader 
is referred to [12] and [13] for an extensive discussion of the eye as a 
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detector. Only general issues relating to the eye and its involvement 
in the simulation will be discussed. 

The eyes are sensitive to positions and velocities. Regan in [56] 
demonstrates that the eyes do not really detect linear accelerations, 
nor do they differentiate velocity to derive acceleration. The eyes are 
most sensitive to velocity in the peripheral vision so the computer 
image generator (CIG) should have a wide field of view, 50 degrees 
or more, to generate adequate velocity cues. Care must be taken at 
these viewing angles to avoid over stimulating the peripheral vision 
if motion is not to be used. These wide angle views tend to induce 
simulator sickness especially when high scene density is present. 

The CIG needs to meet certain requirements. According to Casali, 
[16], the CIG should have proper priority 1 and anti-aliasing hardware 
and software. He also cites other factors contributing to problems in 
visual perception including double imaging or ghosting, and chromi¬ 
nance discontinuities between adjacent CRT’s or screens. When the 
luminance level is too low, it contributes indirectly to perceived dis¬ 
tortion because the pupil will dilate causing an increase in spherical 
aberration. 

Most image systems today are biocular 2 3 rather than binocular z . 
The human vision system is inherently binocular. It is through this 
binocular vision that we perceive depth and distance to objects. A 
biocular CIG system has inherent problems conveying depth and dis¬ 
tance to objects in the display. The eye must be made to focus at 
the proper distance otherwise velocity and distance cueing become 
significant problems. Infinity optics can be used with CRT displays 
to help alleviate this problem, but screen displays have to rely on the 
perspective given to the displayed objects for depth cueing. 

4.1.2.2 The Vestibular System 

The Vestibular system is comprised of the semicircular canals and 
the otoliths. It makes up the non-auditory portion of the ear and is 

1 An assignment to a displayed object which determines whether it is displayed 
in front of or behind other objects. 

2 In a biocular system, both eyes share a common optical axis. 

3 In a binocular system, each eye has its own optical axis. 
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Figure 4.1: The Semicircular Canals 


our primary means to measure linear and angular acceleration. The 
output of these sensors is intimately linked to and correlated with 
the motion our eyes detect. The correlation between vestibular de¬ 
tected motion and visual perceived motion must be maintained in the 
simulator environment for effective simulation. 

The semi circular canals are sensitive to angular acceleration and 
the otoliths detect linear acceleration. There is some sensitivity to 
linear acceleration in the semicircular canals as well as the otoliths to 
angular acceleration. However, the brain is able to resolve this cou¬ 
pling so that there is little problem distinguishing linear from angular 
acceleration. 

Semicircular Canals A typical semicircular canal is shown in Fig¬ 
ure 4.1, taken from Riedel, [60]. There are three semicircular canals 
in each ear arranged in approximately orthogonal axes. Each canal 
is filled with a fluid called endolymph which is similar to water in its 
physical properties. At the base of each canal is a chamber called the 
ampulla. Inside the ampulla is a gelatinous mass called the cupula 
which completely blocks the canal and is attached to nerve fibers 
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through hair like sensory cells. When the head is subjected to an an¬ 
gular acceleration, the endolymph displaces relative to the canal walls 
and in turn displaces the cupula. This in turn changes the afferent 
firing rate of the sensory cells as they bend. The sensory cells in any 
one canal have the same polarization and so all will have a comparable 
change in firing rate for any given displacement, [57,58,60,12,13]. 

The usual model for the semicircular canals is a torsional pendu¬ 
lum. The semicircular canals act like integrating accelerometers and 
integrate angular acceleration to obtain angular rates for normal ve¬ 
locities and durations of head motion. The transfer function for a 
typical model is, [57,58]: 

s~,( \ __ Tl Tg 5 2 _ / , v 

[S} ~ (T l s + 1 )(T S s + 1)(T 0 s + 1) 1 ' ] 

The input is the angular rate and the output is the perceived angular 
rate. The parameters Ti,T a , and 7s are interpreted as time constants 
and s is the Laplace variable. Reid, [57,58], uses the following values; 
10.2 sec., 30.0 sec., and 0.1 sec. respectively. The frequency response 
of this transfer function is given in Figure 4.2. 

The Otoliths The otoliths are separated into a saccule and utricle 
with each ear having one of each. A typical saccule or utricle is shown 
in Figure 4.3, taken from Riedel, [60]. The otolith is a shell filled 
with calcium carbonate crystals, octonia, surrounded by endolymph. 
This system is supported on sensory hair cells similar to those in 
the semi circular canals. When the head is subjected to a linear 
acceleration, the octonia displaces relative to the endolymph changing 
the mass distribution. This change causes the sensory hair cells to 
deflect which, in turn changes the afferent firing rate. 

The utricle is located in a plane that is roughly parallel to the 
horizontal semicircular canal. Its sensory cells are sensitive to motion 
in all directions in that plane. The saccule is located in a plane 
approximately perpendicular to the utricle. However, its sensory cells 
are more sensitive to acceleration perpendicular to the utricular plane, 
[57,58,60,12,13,11] 

The otoliths detect specific force, which is the difference between 
the induced acceleration and the acceleration of gravity. The typical 
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Figure 4.3: The Otolith 

model for the otoliths is a mass balanced accelerometer. Conceptually, 
a mass balanced accelerometer is a device that measures acceleration 
indirectly by the deflection of a beam with a known mass suspended 
at its end. When it is accelerated, an inertial force due to the mass 
causes the cantilever arm to bend. The induced strain of the arm is 
used as an indication of the magnitude of the acceleration. Following 
Reid [57,58], the transfer function model of the otolith is: 

s) = _£_(*f_t!L_ 

[t l s + l)(r s s + 1) 

The parameters t l , r„, and t§ are interpreted as time constants, and 
s is the Laplace variable. Reid, [57,58| assigns the following values; 
5.33 sec., 13.2 sec., and 0.66 sec. respectively. The gain K is given 
the value of 0.4. The frequency response of the otolith model is given 
in Figure 4.4. 

4.1.2.3 Perception Thresholds 

In a simulator, there are certain motions that must be accomplished 
without the subject being aware. For example, if the simulator is to 



Phase (Deg) Magnitude (dB) 


113 




Frequency (Rad/Sec) 


Figure 4.4: Frequency Response of the Otolith Model 
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be tilted to simulate a sustained linear acceleration, the rotation must 
be done at rates that are imperceptible to the subject. A major cause 
of simulator induced sickness is when motion that does not correlate 
with the visual perception is detected by the vestibular system. 

One of the chief problems with assigning a threshold is that there 
is quite a variation in the laboratory results. Examining various ref¬ 
erences such as Berthoz [11], the Handbook of Human Perception and 
Human Performance [12], or Engineering Data Compendium: Human 
Perception and Performance [13] will show this difference. Experimen¬ 
tal information is determined in very controlled circumstances often 
with only one motion axis being excited at any one time. The exper¬ 
iments are also often performed in the dark in soundproofed rooms 
to eliminate any external stimulus. Fortunately, there are other con¬ 
siderations in a simulator which allow us a fair amount of leeway in 
choosing thresholds. In contrast to the above experiments, the simula¬ 
tor environment is actively stimulating many different sensory inputs. 
Although masking of motion in one degree of freedom (DOF) is not 
necessarily masked by motion in another 4 , the high density of infor¬ 
mation from all the subject’s sensors tends to raise all the thresholds. 

In practice, common threshold levels for the semicircular canals are 
in the range of 2.5 to 3.5 degrees/sec and 0.15 to 0.3 meters/sec/sec 
for the otoliths, [11,12,13,57,58]. In the end, knowing the magnitude 
of these thresholds is the starting point. The actual values used are 
best determined by experimentation in a specific simulator given a 
particular experiment. 


4.1.2.4 Proprioceptive and Tactile Sensors 

The proprioceptive sensors are those sensors in the muscles, tendons, 
and joints that respond to stimuli from within the individual. For 
example in the base of the neck there are sensors that help orient the 

4 In recent experiments, Greig, [32] indicates that there may not be a threshold 
per se, but rather that motion detection is better represented by a signal-in-noise 
model and analyzed with signal detection theory. Gaussian receiving operating 
characteristic (ROC) curves are used to define the detector’s performance. The 
detection of motion becomes a statistical question rather than one with absolute 
thresholds. The reader is referred to Greig, [32], for the complete study. 



head. The proprioceptive sensors are normally considered secondary 
to the eyes and vestibular system with regard to motion detection. 

Tactile sensors respond to skin pressure or touch. Again, these are 
secondary in importance to the vestibular system and eyes. 

In certain circumstances, a jet fighter simulator for example, a 
gross motion system is not practical. The accelerations are just too 
large and of too long duration to be simulated with any realism. 
Rather than not having any sense of acceleration, “g-suits” are often 
used. These suits can be pressurized in various ways to stimulate the 
proprioceptive and tactile sensors giving some feeling of acceleration. 


4.1.3 Simulator Sickness 

Simulator induced sickness (SIS) is the target of much concern. A 
simulator rapidly loses its utility if the subjects are forever becoming 
ill during the course of an experiment. Understanding SIS is crucial 
in designing a simulator that will minimize its occurrence. 

SIS shares many common symptoms with motion sickness but is 
more than just a subset of motion sickness. Shared symptoms include 
nausea, pallor, fainting, and a change in respiration and heart rate. 
There are, however, other factors that are exclusively related to SIS. 
For example, SIS can occur when there is no motion at all. Anyone 
who has experienced a very wide field of view motion picture can relate 
to the feelings that this sort of visual stimulation can induce. Visual 
flashbacks are another symptom exclusive to SIS. These flashbacks 
are a concern because they are dangerous. They cause disorientation 
and can occur hours after the simulation has ended. If they happen 
when the subject is for example driving a car, a serious accident could 
occur. 

The causes of SIS are all related to the fact that the simulated en¬ 
vironment is not the same as the real world. Miscues, when vestibular 
and visual motion do not correlate, are a primary cause of SIS. Delays 
in motion onset with regard to the visual motion onset are an example 
of miscueing. The motion of the simulator must be synchronized with 
the visual perceived motion to avoid miscueing. 

Inadequacies in the CIG are suspected as possible causes of SIS 
as well. Casali in [16] lists such things as flicker, off axis viewing, 
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parallax, and image smearing as contributing factors to SIS. 

Sometimes, it’s the high fidelity CIG’s that are more prone to 
inducing SIS. Casali [16] discusses data that demonstrate that a wide 
field of view CIG with a high scene density is more prone to induce 
SIS than a CIG with low scene density when no motion is used. It 
is thought that the reason for this is that the high fidelity system 
stimulates the subject’s vision much more than the low fidelity system 
causing a greater discrepancy between the vestibular system and the 
eyes. If a high scene density CIG is used, motion should almost always 
be used to some degree to help avoid SIS. 

Subjects with a lot of experience in the vehicle being simulated 
are more prone to SIS than those without the same expectations. 
The vestibular system is very adaptable and becomes conditioned to 
expect a certain excitation in familiar circumstances. In a simulator, 
this excitation will often be different. This discrepancy in expected 
and experienced stimulus can result in SIS especially if the exposure 
time to the simulation is extended. 

Masking false motion cues or making up for the lack of motion cues 
is not always possible as Greig, [32], discusses. There has been some 
success, however, found by the army in a recent venture in tank op¬ 
erator training. The army’s SIMNET 5 experiment found that intense 
psychological and aural loading of the subjects in the simulator made 
up nicely for the lack of motion and low grade CIG’s. It’s thought 
that the intensity of the environment occupied the thought processes 
of the participants to the point that the lack of other cues was not 
noticed. 

It is possible that in other circumstances that a highly intense 
emotional environment can mask the lack of correlation in cueing and 
thereby allow more latitude in motion system and CIG design. 

In summary, any difference between the real world and the simu¬ 
lator experience is a possible cause to SIS. It is even more apparent 
when the cues from different sensors do not correlate. The basic phi- 

s SIMNET is a networked tank training environment. There are several simula¬ 
tors linked into one data base. There is no motion and the CIG’s are very basic. 
SIMNET is used to train tank crews in a giant electronic war game. AH the expected 
intense psychological pressure and aural distraction from tank noise, artillery, and 
radio interaction that would be found in a real battle is recreated for the simulation. 
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losophy for avoiding SIS is to recreate all cues as faithfully as possible 
and most importantly, correlate all the cues as well as possible. False 
cues need to be masked or eliminated. If motion contrary to the vi¬ 
sual motion is required, that motion must be at levels imperceptible 
to the subject. 

4.1.4 Designing an Effective Simulator Motion Sys¬ 
tem 

The design considerations of a simulator system are not completely 
obvious nor are they trivial. A complex interaction exists between the 
CIG, motion system, the vehicle’s computer model, and the subject 
in the simulator. This interaction must be considered carefully and 
completely in the design of a system. 

There are always trade-offs in any design effort. A simulator mo¬ 
tion system can be a very large, heavy, and complex system. Payload 
masses of up to 7000 Kg or more are conceivable when one considers 
the mass of a dome, projectors, cab, and other supporting structures. 
When this mass is supported on a motion system, large moments of 
inertia can result. These large masses and moments of inertia cre¬ 
ate problems for sizing the hydraulics and compensation in the servo 
controller. Careful thought and design has to be used in order to 
minimize the cost and physical size of the simulator and hydraulic 
systems. 

Another trade-off is related to the motion required for the ve¬ 
hicle being simulated. The platform performance envelope must be 
matched carefully to the performance envelope of the vehicle being 
modeled. An aircraft has different predominant motions than a land 
based vehicle. The maneuvers are different and the control is different. 

A motion system that is optimal for an aircraft is most certainly sub- 
optimal for an automobile. A careful evaluation of the performance 
and characteristics of a vehicle must be done prior to the design a 
motion system to simulate that vehicle. 

Even when the simulator motion system is optimal, it will still 
not be able to exactly replicate all the motions of the vehicle. The 
output of the vehicle dynamics needs to be processed so that the 
motion system is commanded in such a way as to make the best use 
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of the available performance at the same time giving the proper cues 
to the subject in the cab. A washout filter is used for this purpose. 
A fundamental piece of the washout is a high pass filter. The high 
pass function removes the low frequency components of the motion 
commands that are the most likely to cause performance limits to be 
exceeded. Motion commands that have a frequency content typically 
of 2 Hz or lower are those that cause the most problems. 

Washout filters deserve special attention because of their key role 
in the effectiveness of a simulator. 

4.1.4.1 Washout Filters 

A washout filter is more than just a digital filter. It provides rate 
and acceleration limiting when needed to keep certain motions below 
thresholds, as well as providing the means of simulating those motions 
that are beyond the capabilities of the motion system. Washout filters 
play a critical role in how effective the simulator’s motion system is. 

Tilt Coordination The hardest motion to simulate is a sustained 
acceleration. Eventually, any simulator motion base will run out of 
displacement capability and sustained accelerations quickly push a 
simulator to its displacement limits. The most common method used 
to simulate sustained accelerations is with tilt coordination. Sustained 
lateral and longitudinal accelerations are just not feasible using most 
standard simulator motion bases. An obvious solution is to tilt the 
simulator and let gravity make up for the lack of performance. This 
is the basic principal of tilt coordination. There is one significant 
difficulty with this, if the simulator is tilted too fast, the subject will 
be miscued and could succumb to SIS. 

Tilt coordination has had good success in flight simulators because 
of the very nature of the maneuvers performed in an aircraft. In gen¬ 
eral, only large relatively slow response aircraft use motion systems. 
This slow response allows time for tilt coordination whereas in an 
automobile, the response is much faster and the maneuvers are more 
sudden. The onset and termination of accelerations are very fast so 
there is just not time for proper tilt coordination. 

Imagine coming to a sudden stop. The simulator will typically run 
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out of displacement and tilting the simulator is a natural solution. 
However, when the brake is released, the acceleration is terminated 
almost instantaneously. If the simulator remains tilted, the subject 
still feels that he is stopping. If the tilt coordination is terminated 
quickly, the subject will perceive the rotation and a miscue will result. 
Hahn and Kalb in [27] found that tilt coordination was a problem in 
their experiments in the Daimler Benz Simulator. They suggest that 
information from the visual data base could be used to anticipate the 
possibility of sustained lateral and longitudinal accelerations. These 
data could be used to intelligently make use of the simulator’s perfor¬ 
mance envelope. Normally though, tilt coordination is not useful in 
a driving simulator. 

The correlation of the tilt coordination is achieved through a low 
pass filter. The low frequency and constant linear accelerations need 
to be coupled to the appropriate rotational DOF. The lateral acceler¬ 
ation is coupled to the roll and the longitudinal acceleration is coupled 
to the pitch. There is no coupling between the vertical acceleration 
and the yaw. The low pass filter is designed such that the convolution 
of its transfer function and the associated high pass filter’s transfer 
function is approximately unity. In this way, the frequency content of 
the resulting motion command for a linear acceleration is representa¬ 
tive of the frequency content of the original signal. Rate limiting of 
the tilt coordination signal is used in order to avoid the perception of 
a miscue. 

Another method to sustain linear accelerations is conceptualized 
in [36]. The proposed motion system has a continuous yaw capabil¬ 
ity. This would allow the simulator to be run in a centrifuge mode 
so that constant lateral, or in certain instances longitudinal, acceler¬ 
ations could be maintained indefinitely. The washout filter can de¬ 
signed to implement this sort of mode instead of tilt coordination. 
The low passed filtered specific force can be used to determine certain 
parameters regarding the centrifuge mode such as the radius and rate 
of rotation. 

Scaling and Limiting The washout will usually include some means 
of scaling the output of the vehicle dynamics. It is common to reduce 
the magnitude of the input to the washout by as much as half in 



120 


CHAPTER 4. PHYSIOLOGY FOR SIMULATION 


flight simulators. Attenuating the input by more than this will tend 
to reduce the effectiveness of the resulting motion, [64]. The effect 
of scaling is to effectively increase the performance envelope of the 
motion system by the inverse of the scale factor. The risk is that the 
subject may perceive that the accelerations are not of the magnitude 
that would be expected. 

There are usually some sort of limit switches in hydraulic systems 
that will shut down the system if limits are exceeded. The washout 
often contains software limiting [57,58] that will avoid these limits so 
that panic shutdowns are avoided. 

Two approaches to washout will be further investigated. The first 
is based on classical linear filter theory and is used to provide the 
basic frame work for conceptualizing a washout filter. The second 
is a filter based on an adaptive concept that overcomes some of the 
weaknesses of the classical linear filter. 


Linear Filter Theory The development of a standard linear washout 
filter follows conventional digital filter design. In practice, a second 
or third order high pass filter is used. 

The basic data flow for the linear filter is shown in Figure 4.5 
following, [57,58,10,47,60]. The input to the filter is the specific force 
and the angular rate. The output is the high pass filtered specific 
force and a combination of the high pass filtered angular rate and low 
pass filtered specific force. The transfer function representation of a 
typical high pass filter for the translational DOF is, [57,58]: 


F(s) = 


(s’ + 2 ftp WAp S + ujp)(s + U>Ap) 


(4.7) 


Where Uhp is the cutoff frequency and is the damping factor. 
The complementary low pass filter is given by: 


r(*) = 




IP 


S 2 + 2 G P u tp s + wf p 


(4.8) 


Where a, = ft, and w lp = 2 u kp 
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The resulting equations are commonly converted to the discrete 
domain by use of the bilinear z transform 6 or by integrating the time 
domain equations numerically. 


Adaptive Washout Filter The linear filter works well but can 
generate spurious false cues that can be very disconcerting to the 
subject in the simulator. A conceptual drawing showing the false 
cue is shown in Figure 4.6. Ariel and Sivan, [3], and Reid and Nahon, 
[57,58], present an adaptive filter structure that effectively reduces the 
false cue to imperceptible magnitudes. Parrish, (55,60) has a scheme as 
well but Ariel and Sivan show that it can distort the motion perception 
for some types of motion commands. The filter presented here follows 
the development by Reid and Nahon, [57,58]. 

The intent of an adaptive washout filter is to adapt the severity of 
the filter to the current state of the simulator. If there remains sig¬ 
nificant motion capability in the simulator, the filter need not restrict 
the motion commands at all. When limited capability remains, the 
filter needs to drastically attenuate all motion commands. 

The design of an adaptive filter starts with defining the differential 
equation that relates input and output. The general form of the 
adaptive filter equation for translation is taken as [57,58,3]: 

iJL = Pl a - kiS - (4.10) 

Where: 


• S = Displacement vector from the neutral position to the cur¬ 
rent simulator position 


• a = Acceleration of the vehicle 
G The bilinear z transform is: 


- 2 ( 1 ~ z ~ 1 ) 
5 rii + z- 1 ) 


(4.9) 


Where T is the sample period. The advantage of the bilinear z transform is that it 
maps the left hand side of the s-plane completely into the unit circle of the z-plane. 
The disadvantage is that there is a frequency warping effect that must be taken into 
consideration, (53). 
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Figure 4.6: False cue generated with a step input to the high pass 
filter. 
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• pi = the adaptive parameter 

• k u k 2 = constant weighting factors 

For rotation [57,58,3]: 

^ = t{p 2 a) +Ps^ (4.11) 

Where: 

• V, = angular velocity of the simulator 

• u> c = angular velocity of the vehicle 

• P 2 iPz = adaptive parameters 

• a = acceleration of the vehicle 

• p 2 a = A term that represents tilt coordination. 

• t Refers to a non-linear limiting function that is applied to this 
particular term. It performs no limiting until the term reaches 
a pre-determined value then the limiting function clamps the 
value of this term at this limit. This clamping is maintained 
until the value falls below the limit. 

A cost function is designed next. In this case, the cost function in¬ 
cludes terms that define errors that should be minimized, and penalty 
terms that help increase the filtering effect in certain regions. The 
method of steepest descent is used to define the rate of change of the 
adaptive parameters. The cost function is differentiated with respect 
to the adaptive parameters and the resulting equations are used with 
the original differential equations to form a set of equations that are 
then solved numerically. 

Differentiating the cost function has the effect of optimizing the 
filter equation for a given state of the input. The adaptive parameters 
are continually being updated so the filter is always tending towards 
an optimal design. 

Here is a simplified example of the cost function for a translational 

DOF: 
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J, = r[(<* — 5) 2 +(^cor — ^»m)* + (5*+S 2 +^J im H-^ m )-t-(pi-pi [ ,) I +(pi—pj tf ) 2 + (p3—pg#) 2 ] 

(4.12) 

Where: 

• a — 5 = the error between the acceleration of the car and the 
simulator. 

• <t>c or — 4>iim = the error between the roll velocity of the car and 
the simulator. 

• S 2 + S 2 + 4>\ irn + = A penalty weighting factor based on 

the linear velocity, linear position, angular velocity, and angular 
position of the simulator. 

• Pn ~ Pno — The difference between the initial value of the adap¬ 
tive parameter and the current value. 

4.1.4.2 Summary 

In summary, the design of an effective simulator must weigh the re¬ 
quired performance against the physical limitations of the technology 
at hand. Where the motion system lacks performance, proper use of 
washout can compensate to some degree. It’s also possible to make use 
of other stimuli such as intense psychological loading, to compensate 
for limits in other systems. 

The overall system must be matched to the expected maneuvers 
as well as the general vehicle performance envelope. Simply choos¬ 
ing the highest performance subsystems and connecting them is not 
enough. A balance must exist between the motion, visuals, and sound 
to produce the most ideal simulated environment. 

4.1.5 Control Force Loading 

An effective simulated environment does more than recreate visual, 
aural, and gross motion stimuli. There are other ways in which we 
interact with our environment when driving or flying. In particular, 
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we receive feedback through the controls that we use to maneuver a 
vehicle. The reaction forces and torques we are subjected to need to 
be recreated in the simulator environment in a realistic way. 

There also has to be a link from the driver’s controls to the input 
of the vehicle dynamics. 

4.1.5.1 Vehicle Commands and Controls 

The vehicle dynamics model needs to know what the driver is doing. 
Data acquisition is performed on the steering, brake, throttle, gear 
selection and any other controls that are available to maneuver the 
vehicle. Acquiring the command data is rather straight forward. Basic 
instrumentation is used and adapted to the needs of the particular 
model. These data are then made available to the vehicle dynamics 
at regular intervals. 

The driver of a vehicle expects to have a certain reaction force on 
each of the controls. Sometimes this reaction is nothing more than 
a spring such as a throttle pedal. At other times, the force is linked 
to the dynamics of the vehicle such as for the steering. The correct 
forces need to be determined and then applied to the various controls. 

Determining the control force for controls such as a throttle or 
clutch is straight forward. These systems are all spring loaded and 
adequate feedback is obtained from simply loading these controls with 
the appropriate springs and dampers. 

The steering is a harder problem. The driver will always have a 
hand on the wheel and uses it to derive a lot of information about ve¬ 
hicle handling and road feel. This makes proper loading of the steering 
wheel a very important issue to proper and effective simulation. 

The driver reacts forces and torques from the tires through the 
steering system. These forces and torques include an aligning torque 
from the wheel caster, lateral forces from cornering, and forces and 
torques from the tire-road interaction. This means that the vehicle 
dynamics needs to calculate these reactions as part of the regular cycle 
and communicate the resulting acceleration of the steering wheel to a 
control system for the steering. 

In the particular dynamics formulation that we use, [22], the steer¬ 
ing rack is modeled as one of the bodies. The acceleration, velocity 
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and displacement of the rack are then available for use at the end of 
each integration. The acceleration can be used as the setpoint for a 
servo control system connected to the steering wheel. The reaction 
to this acceleration command is then provided by the driver just as it 
would be in a real vehicle. 

The transmission can be straight forward especially if it is an au¬ 
tomatic transmission. In this case, no real feedback is present other 
than the shift lever position. This can be handled by the original 
shifting mechanism or by detente springs. 

A manual transmission presents different problems. There is more 
feedback particularly if it is a floor mounted shift lever that enters the 
transmission directly. A fairly complicated control system is needed 
unless the real transmission can be used. Even then, the transmission 
must be driven at the input and output so that gear grinding and the 
proper feel when shifting is achieved. Fortunately, the transmission 
is used intermittently and the control force feedback is not critical 
to the simulation. This allows a lot of latitude in the transmission 
feedback system design and implementation. 

The gages need to be driven as well. Parameters such as the speed, 
and rpm are available from the dynamics and drive train models. The 
oil pressure, water temperature, battery charging, and any other gages 
can normally be set to a nominal value and left alone. 

4.1.6 The Total System 

Integrating the simulator sub-systems properly presents a great chal¬ 
lenge. Each of the separate parts interacts intimately with the whole. 
It is not enough to design each sub-system to perform to its maximum 
capability. Rather, it is better to design each component so that it 
functions with comparable effectiveness as all the rest. A bottom 
up design approach where pieces are put together without regard to 
overall system performance requirements does not work. The overall 
system performance must be considered and broken down into pro¬ 
gressively smaller systems. This top down approach will insure that 
the final system will perform as wanted. The example of using a flight 
simulator motion base for a driving simulator is a good example. A 
flight simulator motion system has motion capabilities in keeping with 
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the motion of typical aircraft. The same motion base will be wholly 
inadequate for a land based vehicle. Its capabilities may exceed the 
needed motion in some DOF’s for an automobile but will fail miser¬ 
ably for others. 

An effective simulator system will be balanced mix of sensory stim¬ 
ulating devices. A system with modest capabilities in vision, motion, 
sound, and control loading will in general produce a better simulated 
environment than one that has a very expensive motion system at the 
expense of the CIG, sound and control loading. 

4.1.7 The Future 

The future of vehicle simulation is leading towards higher performance 
CIG’s and vehicle dynamics. As hardware performance improves, it 
will open the door to CIG’s that have higher scene detail, faster frame 
rates, and shorter pipeline delays. The vehicle dynamics will be able 
to take advantage of the faster throughput to include more detail 
in the modelling. For example, the addition of non-linear bushings, 
flexible bodies, and more complicated analytical tire models are but 
a few of the areas of current research. 

The future holds great promise for the fidelity of real time sim¬ 
ulation. New hardware, algorithms carefully designed to make use 
of parallel processing and new computing power, and more efficient 
mechanical systems will provide the means to achieve greater realism. 
The utility of simulators for training and engineering can do nothing 
but increase as performance goes up and cost goes down. 



Chapter 5 
References 


5.1 References on Dynamics 

Because the current vehicle models decomposes the car, or object or 
whatever into a set of linked rigid bodies, we must remember how the 
dynamics and kinematics of a single rigid body is analyzed. The best 
single reference on most topics in theoretical mechanics is Goldstein’s 
book [30). Beware that many misprints exist in the text and exercises. 


An important reference for the foundations of the virtual work prin¬ 
ciple and recursive dynamics is Haug’s excellent soon-to-be-published 
(?) text (37). It is a very readable treatment on a less Olympian level 
than Goldstein, and even if less advanced (no treatment of Hamil¬ 
tonians equations), it is more directly useful in understanding the 
tree-based recursive vehicle dynamics formulation. Haug, as well as 
his students, use the notation and terminology found in his text. 

A good short treatment of dynamics 1 is Landau’s book [45] of some 
hundred pages, each of which could be written in stone. Finally, 
among the mechanics texts, we must mention Lanczos’ book [44] on 
variational principles in mechanics. It has has few formulas, but much 
insight, and it reads like a novel. In addition to these gems, there are 
many engineering and physics mechanics texts available to match any 

1 Based almost exclusively on Lagrangian methods. 
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taste. 

Because at times we’ve mentioned the underlying geometric and topo¬ 
logical nature of mechanics, we’ve listed a few decent geometry refer¬ 
ences for the non-mathematician. To begin, an excellent text devoted 
to differential geometry for the scientist is the recent work by Burke 
[14). The book is dedicated to “a// those who, like me, have wondered 
how in hell you can change q without changing q " (an irresistible 
invitation). Reference [19] is a very complete treatment of modern 
methods in mathematical physics and includes an excellent section on 
differential geometry and manifolds. Another text on things differ¬ 
ential geometric 2 is the famous older book by Synge and Schild [67]. 
We recommend avoiding the large tome by Abraham and Marsden 
[l], except possibly to look up specific topics once you are familiar 
with the heavy machinery surrounding all the physics. 


2 But from the classical “tensors as indexed objects” approach, thereby avoiding 
the more modern global methods. 
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5.2 References on Perception 

The concept of washout is probably new to many readers. A good 
place to start is with Reid and Nahon, [57,58]. Thewe reports pro¬ 
vides an excellent introductino to washout and human perception 
modelling. Riedel’s, [60], work is similar but has differences in the 
exact formulation of the washout and perception models. Both ref¬ 
erences provide an excellent groundwork for understanding washout 
and related issues to simulator motion concerns. 

Much work has been done in recent years regarding the physiology 
of the human perception systems. The two works, [12,13], edited by 
Boff are extremely detailed volumes on the physiology and operation 
of human perception. Regarding the actual modelling of vestibular 
systems, Riedel, [60], or Reid and Nahon, [57,58], are starting points 
for basic modelling issues. Greig, [32], has some new ideas, partic¬ 
ularly regarding threshold levels, that are worth study by interested 
readers. 
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